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CHAPTER  I 


INTRODUCTION 


1.1  Problem  Description 

The  restoration  of  astronomical  images  to  permit 
greater  resolution  is  a problem  of  long-standing  interest. 
Image  resolution  is  limited  both  by  atmospheric  turbulence 
and  by  the  diffraction  effects  caused  by  imaging  Kith  a 
finite  aperture,  although  the  former  is  more  serious  for 
ground-based  telescopes. 

Atmospheric  turbulence  can  be  thought  of  as  causing  a 
random  perturbation  of  the  refractive  index  of  air  in  the 
optical  path.  This  perturbation  varies  in  both  time  ard 
space,  and  causes  the  phase  of  an  incident  light  wavefront 
to  be  distorted.  These  distorted  waves  are  not  brought 
into  sharp  focus  by  a conventional  imaging  system,  but  are 
instead  spread  over  a large  area  in  the  image  plane.  For 
the  200  inch  Hale  telescope,  the  diameter  of  the 
diffraction-limited  Airy  disk  in  visible  light  is  about  .06 
arc-seconds.  However,  fo>'  observations  made  through  the 
atmosphere  the  minimum  signal  diameter  in  the  image  plane 
is  about  2 arc-seconds,  which  corresponds  to  a 40-fold 
decrease  in  resolution. 

For  many  years  this  lost  resolution  was  considered 


?. 

i r reco  verab  1 e . However,  in  the  early  1970s  i_abeyrie  showed 
that  di f f r ac t i on- 1 i m i ted  information  could  be  obtained  by 
ground-based  telescopes.  His  approach  was  to  take  a series 
of  very  short  exposure  photographs,  each  of  which 
essentially  'freezes'  the  motion  of  the  atmosphere. 
Because  of  the  random  constructive  and  destructive 
interference  of  the  distorted  light  waves,  the  resulting 
photographs  have  a speckled  appearance.  If  the  Fourier 
transform  of  each  photograph  is  taken,  the  phase  spectrum 
is  seen  to  vary  randomly  from  photograph  to  photograph. 
When  these  short-exposure  photographs  are  averaged  directly 
(which  is  equivalent  to  long-exposure  imaging),  the  random 
phase  spectrum  variations  cause  cancellation  of  the  high 
spatial  frequencies.  The  result  is  a 1 ow- reso 1 uti on  imaqe 
typical  of  conventional  astronomy.  To  avoid  this  result, 
Labeyrie  ignored  the  phase  component  and  averaged  only  the 
squared  modulus  of  the  Fourier  transforms.  The  resulting 
estimate  is  a diffraction-limited  estimate  of  the  power 
spectral  density  of  the  image,  or  equivalently  the  image 
autocorrelation.  Because  no  estimate  of  the  phase  spectrum 
is  obtained,  however,  an  image  of  an  arbitrarily  shaped 
object  cannot  be  formed. 

To  supply  the  phase  spectrum  estimate,  Knox  and 
Thompson  proposed  a new  technique  for  processing  the  same 
speckle  photographs.  They  multiply  the  Fourier  transform 
of  each  photograph  by  a shifted  copy  of  itself  and  average 
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this  product  over  the  series  of  photographs.  This  averaged 
product  estimates  phase  differences  between  points  in  the 
spatial  frequency  domain,  which  can  then  be  integrated  to 
form  a phase  estimate.  Knox  indicates  that  the  estimator 
is  very  sensitive  to  sensor  noise  in  the  observations. 

Although  encouraging  results  have  been  obtained  in 

simulations,  no  successful  reconstructions  of  real 

astronomical  objects  have  been  obtained  to  date.  y j 

As  mentioned,  diffraction  effects  also  limit  the  p;, 

» ij 

possible  image  resolution  by  imposing  an  absolute  spatial 
frequency  cutoff  on  the  signal  spectrum,  called  the 
diffraction  limit.  This  limit  is  a function  of  the 
wavelength  of  the  observed  light  and  the  diameter  of  the 
telescope  aperture.  In  the  absence  of  telescope 
aberrations  the  optical  transfer  function  (OTF)  of  the 
telescope  is  essentially  a triangular  function  that 
low-pass  filters  an  input  signal.  A number  of  techniques 
have  been  proposed  to  extrapolate  the  measured  spatial 
frequency  spectrum,  the  most  successful  being  the  iterative 
algorithm  of  Gerchberg.  This  algorithm  alternatively 
imposes  constraints  in  both  the  spatial  and  frequency 
domains  to  perform  the  extrapolation.  As  Gerchberg  himself 
pointed  out,  the  algorithm  tends  to  corrupt  the 
extrapolated  spectrum  at  each  iteration  by  adding  to  the 
extrapolated  spectrum  a portion  of  the  distortion  energy  in 
the  measured  spectrum.  Thus,  extrapolations  of  noisy  data 
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obtained  by  this  algorithm  may  he  suspect.  Despite  this 
caveat,  the  algorithm  has  recently  been  used  to  aid  In  a 
digital  recons t rue 1 1 on  of  Retelgeuse. 

These  two  problems,  that  Is,  accurately  estimating 
phase  spectra  and  extrapolating  low-frequency  spectra,  are 
common  to  a great  many  deconvolution  problems.  There  Is  to 
date  no  entirely  satisfactory  solution  to  this 
deconvolution  problem  in  the  general  case.  The  fundamental 
difficulties  are  those  of  Inversion  theory;  namely,  the 
truncation  of  Infinite  impulse  responses,  smoothing,  and 
noise  amplification.  In  the  blind  deconvolution  problem,  a 
spectral  estimation  problem  must  also  be  solved.  Phase 
estimation  has  traditionally  been  more  difficult  than 
magnitude  estimation,  so  much  so  that  researchers  often 
ignore  it  when  dealing  with  speech  and  music  signals,  where 
it  is  often  not  perceptually  significant.  However,  for 
Images,  such  a cavalier  approach  leads  to  poor  results. 
Thus  the  astronomical  imaging  problem  presents  an  important 
special  case  of  the  deconvolution  problem  because  of  the 
poor  SNR  and  the  severity  of  the  phase  distortion. 

1.2  Contribution  of  this  Research 

This  research  describes  our  efforts  to  Improve  the 
resolution  of  astronomical  Images,  first  by  modifying  the 
Knox-Thompson  estimator  to  provide  an  adequate  phase 
spectrum  estimate  out  to  the  telescope  diffraction-limit, 
and  second  by  estimating  the  spatial  frequency  spectrum 
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beyond  the  d i f f r ac ti on- 1 im 1 1 . The  techniques  described 
here  have  been  strikingly  successful,  and  we  anticipate 
that  they  will  be  applied  to  a variety  of  deconvolution 
problems. 

Our  most  important  contribution  to  phase  estimation  is 
the  application  of  a fast  direct  algorithm  to  perform  a 

« 

least- squares  int eg  ration  of  the  estimated  phase  ( 

' 

differences.  The  algorithm  is  stable,  accurate,  ^ 

k 

memory-efficient,  and  about  an  order  of  magnitude  faster 

M 

than  the  best  iterative  schemes.  The  resulting  estimates 
are  more  accurate  than  those  obtained  by  integrating  the  f* 

phase  differences  along  several  different  paths  of  i 

integration  and  then  averaging  these  discrete  line 
integrals.  We  show  empirically  that  even  better  estimates  i 

can  be  obtained  by  refining  the  least-squares  phase  with  a 
Gerchberg- Sax  ton  algorithm.  The  practical  necessity  of 
windowing  the  speckle  Images  is  also  demonstrated. 

The  other  principal  contribution  Is  the  development  of 
a new  super- resol uti on  algorithm,  based  on  a non-linear 
modification  of  the  Gerchberg  algorithm.  This  new  approach 
is  more  accurate  in  the  presence  of  noise  and  demonstrates 
superior  convergence  properties  when  compared  to  the 
original  al gori thm . 

These  techniques  have  been  used  to  restore  images  of 
the  red  super-giant  star  Betelgeuse  and  the  binary  star 
Capella.  These  are  the  first  successful  reconstructions  to 
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explicitly  estimate  the  object  phase  spectrum.  The 
reconstructions  are  consistent  with  known  a s trophy s i c al 
data  and  appear  to  be  of  high  quality. 

This  dissertation  is  arranged  as  follows.  The 
previous  research  on  the  problems  of  speckle  imaging,  phase 
estimation,  and  super- resol uti on  are  reviewed  in  chapter 
II.  Chapter  III  describes  in  detail  the  1 east- squares 
phase  integration  algorithm,  and  compares  it  to  the  line 
integral  estimator  via  simulations  for  various  SNR.  The 
use  of  various  modifications,  including  phase  refinement, 
is  also  explored.  Chapter  IV  is  a similar  treatment  of  the 
new  super- reso 1 uti on  algorithm.  Chapter  V describes  the 
experimental  technique  and  the  results  obtained  using  real 
images.  A brief  summary  of  the  results  is  reviewed  in 
chapter  VI,  including  some  remaining  areas  of  research. 


CHAPTER  II 


♦ 


LITERATURE  REVIEW 

II. 1 Speckle  Interferometry  and  Some  Extensions 

The  difficulty  in  imaging  through  the  turbulent 
atmosphere  has  long  been  recognized,  and  astronomers  have 
resorted  to  the  use  of  interferometers  to  obtain  accurate 
quantitative  data  on  the  angular  dimensions  and  separations 
of  stellar  objects.  The  classical  astronomical 
interferometers,  the  Michelson  stellar  interferometer  and 
the  intensity  interferometer,  invented  by  Hanbury  Brown  and 
Twiss,  are  inherently  one- dimen s i onal  devices,  however. 
They  are  therefore  most  useful  when  a one- dimensional  scan 
sufficiently  characteri zes  the  object  of  interest.  For 
extended  objects  pictorial  information  is  often  desired, 
and  one-dimensional  techniques  do  not  suffice. 

The  development  of  speckle  interferometry  by  Labeyrie 
[II  and  his  successful  use  of  the  procedure  [2]  to  measure 
object  autocorrelations  and  hence  to  infer  angular 
dimensions  represent  a step  towards  achieving  full  aperture 
reconstructions.  The  method,  which  averages  only  the 
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squared  modulus  of  the  Fourier  transform  of  the  image,  is 
successful  precisely  because  it  ignores  the  phase.  This 
prevents  the  cancellation  of  the  high  spatial  frequency 


energy,  but  it  also  precludes  the  reconstruction  of  an 
image.  As  such  it  represents  only  a partial  solution  to 
the  problem  of  d i f f r ac ti on- 1 im i ted  imaging. 

Since  the  publication  of  Labeyrie's  work  a number  of 
researchers  have  proposed  extensions  to  the  technique  to 
permit  true  imaging.  Most  of  these  proposals  avoid  the 
problem  of  phase  estimation  by  working  directly  with  either 
the  speckles  themselves  or  the  estimated  object 
autocorrelation.  The  earliest  suggestions  (those  of  Bates, 
Gough  and  Napier  [3]  and  Liu  and  Lohmann  [41)  required  the 
use  of  an  unresolvable  point  star  somewhere  in  the  same 
isoplanatic  patch.  ( I sopl anic i ty  is  synonymous  with 
shift-invariance).  This  is  a rather  unsatisfactory 
criterion,  since  it  is  often  not  possible  in  practice. 
Weigelt  [5]  simulated  a variation  of  the  Liu  and  Lohmann 
technique,  wherein  he  varied  the  brightness  of  the 
reference  star  to  determine  which  of  the  apparent  stars  in 
the  image  autocorrelation  were  genuine  and  which  were 
artifacts.  The  artifacts  were  then  removed 
photographically  during  the  final  copying  process.  The 
relative  brightness  of  the  reference  star  was  varied  by  the 
use  of  a neutral  density  filter  placed  in  front  of  each 
speckle  image.  In  fact,  the  filter  was  used  to  darken  the 
speckles  which  are  not  produced  by  the  reference  star. 
Weigelt  reported  good  results  in  his  simulation,  even  when 
the  object  field  was  so  large  that  it  no  longer  lay 


entirely  within  an  isoplanatic  patch. 

More  recently,  some  interesting  results  have  been 
obtained  from  real  speckle  images  by  an  entirely  different 
approach.  Lynds,  Worden  and  Harvey  [6]  describe  a 
technique  in  which  it  is  assumed  that  the  speckles  are 
formed  by  the  convolution  of  the  true  image  with  a 
point-spread  function  consisting  of  an  array  of  Dirac  delta 
functions,  each  of  which  is  located  at  the  center  of  mass 
of  a speckle.  The  problem  is  then  one  of  deconvolution. 
An  editing  procedure  selects  the  brightest  and  clearest 
speckles  so  as  to  most  accurately  estimate  the  location  of 
the  delta  functions.  The  deconvolved  result  is  then 
averaged  over  a number  of  images  to  improve  the  SNR.  They 
carried  out  a demonstration  of  this  procedure  on  speckle 
photographs  of  the  red  supergiant  Alpha  Orionis 
( Betel geuse)  , and  indeed  found  some  evidence  of  fine 
surface  structure.  There  were  internal  evidences  that  the 
technique  did  not  achieve  full  diffraction-limited 
resolution,  however,  despite  the  increased  image  detail. 

To  improve  the  resolution  of  the  Lynds  et  al  . result, 
a further  refinement  was  implemented  by  McDonnell  and  Pates 
[7],  They  accept  the  result  as  a preliminary  estimate,  but 
assume  that  an  array  of  simple  delta  functions  is  an 
oversimplification,  and  that  the  resulting  image  is  still 
blurred,  although  to  a lesser  extent,  by  a residual 
point-spread  function  (PSF).  The  estimation  of  this 
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residual  PSF  is  rather  involved.  The  Lynds  et  al  . result 
is  assumed  to  have  as  the  zeros  of  its  Fourier  transform 
the  union  of  the  zeros  of  the  Fourier  transforms  of  the 
object  and  the  residual  PSF.  To  separately  identify  these 
zeros  an  intensity  model  for  the  object  is  assumed,  that  of 
a simple  first  order  limb-darkened  star.  The  identity  of 
the  image  zeros  is  determined  by  comparison  with  reasonable 
values  predicted  by  the  model.  The  zeros  of  the  residual 
PSF  are  thus  identified,  and  the  image  is  Weiner  filtered 
by  the  inverse  residual  PSF.  The  result  is  an  image  whose 
size  is  very  near  the  telescope  diffraction  limit,  and 
certain  astrophysical  features  are  visible. 

This  restoration  technique  has  the  advantage  that  no 
observation  of  a reference  star  is  required.  The  Lynds  et 
al . assumption  only  holds  for  small  bounded  objects  whose 
angular  size  is  about  the  same  as  the  telescope  diffraction 
limit,  which  restricts  the  usefulness  of  the  technique 
somewhat.  The  McDonnell  and  Bates  refinement  would  be 
difficult  to  extend  to  objects  with  non-simple  intensity 
distributions,  because  a more  complex  set  of  zeros  might 
defy  an  easy  segregation  of  PSF  and  object  transform  zeros 
into  separate  sets.  In  addition,  the  need  to  assume  an 
intensity  distribution  model  begs  the  question  since  the 
intensity  distribution  is  precisely  what  we  are  trying  to 
find. 

A related  processing  technique  which  seeks  to  remove 
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the  necessity  of  point  star  observation  is  suggested  by 
Welter  and  Worden  [ 8 J . They  show  that  the  difference 
between  the  average  autocorrelation  of  instantaneous 
speckle  images  and  the  crosscorrel ation  of  the  same  images 
is  proportional  to  the  image  autocorrelation  undegraded  by 
atmospheric  turbulence.  The  approximation  is  most  accurate 
if  the  observed  object  is  small  compared  to  the  seeing 
spot.  How  this  technique  may  apply  to  speckle  imaging  is 
not  yet  known. 

Sherman  [9]  describes  a technique  whose  application  is 
unhindered  by  restrictions  on  the  location  of  a point  star 
reference  or  the  object  size,  so  long  as  the  image  lies 
within  an  isoplanatic  patch.  He  proposes  the  calculation 
of  the  entire  four-dimensional  covariance  function  of  the 
Fourier  transform  of  the  image.  The  object  transform  is 
then  estimated  as  the  complex  eigenmatrix  associated  with 
the  largest  eigenvalue  of  the  covariance.  He  obtains  very 
good  results  in  a one-dimensional  simulation,  and  the 
technique  is  claimed  to  be  much  more  immune  to  measurement 
noise  than  the  Knox-Thompson  method  discussed  in  the  next 
section.  The  amount  of  computation  needed  to  implement 
completely  this  technique  would  make  it  impractical  for  all 
but  the  most  powerful  computing  facilities.  In  T10] 
Sherman  describes  a modification  which  greatly  reduces 
computation.  At  the  present  time,  his  method  appears  to  be 
the  most  practical  alternative  to  explicit  phase  estimation 
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for  a general  image. 

II. 2 Explicit  Phase  Estimation 

Phase  estimation  has  remained  an  extraordinarily 
difficult  problem  in  image  restoration.  Recent  work  by 
Cole  [11]  demonstrates  the  rather  poor  quality  restorations 
which  result  when  real  world  blurs  are  assumed  to  have  zero 
phase.  Cannon  [12]  showed  some  improvement  by  the  simple 
expedient  of  considering  the  phase  to  be  either  0 or  pi, 
that  is,  phase  reversals.  However,  the  first  attempt  to 
actually  estimate  the  image  phase  appears  to  have  been  done 
by  McGlamery  [13].  He  averaged  separately  the  magnitude 
and  phase  spectra  of  each  of  the  speckle  images.  He 
concluded  that  as  the  variance  of  the  phase  estimate 
approaches  pi/3,  the  phase  information  provided  by  this 
estimator  becomes  vanishingly  small. 

The  technique  of  central  concern  to  this  research  is 
that  of  Knox  and  Thompson  [14].  They  proposed  estimating 
the  image  phase  by  computing  the  one-shift  term  of  the 
autocorrelation  of  the  Fourier  transform  of  each  speckle 
image.  This  term  preserves  the  phase  information  in  the 
form  of  phase  differences  in  the  x-  and  y-directions 
between  points  in  the  Fourier  transform.  These  phase 
differences  must  then  be  integrated  to  form  a phase 
estimate.  In  his  Ph.O.  dissertation  Knox  [15]  further 
showed  by  both  analytical  means  and  simulations  the 
relative  insensitivity  of  the  technique  to  fixed  telescope 
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aberrations.  This  is  fortunate,  since  it  may  permit  the 
use  of  this  technique  in  conjunction  with  large  multiple 
mirror  tele sc  opes,  which  often  have  substantial  optical 
aberrations.  A similar  extension  of  speckle  interferometry 
has  already  been  worked  out  by  Roddier  [16].  Knox  also 
found  that,  due  to  the  narrowness  of  the  random 
instantaneous  optical  transfer  function  (OTF),  the 
normalized  SNR  of  the  terms  in  the  autocorrelation 
decreases  rapidly  as  the  shifted  distance  increases,  so 
that  at  best  only  the  first  few  terms  contribute 
significant  phase  information.  Knox  suggested  that  the 
phase  be  estimated  by  averaging  the  phase  estimates  derived 
from  summing  the  phase  differences  over  two  independent 
paths.  His  one-dimensional  simulations  showed  that  very 
good  results  could  be  obtained.  They  also  showed  that  the 
final  reconstruc ti ons  were  quite  sensitive  to  errors  in  the 
phase  spectrum  estimate. 

At  Itek  Corporation  a two-dimensional  simulation  of 
the  Knox-Thompson  (K-Tl  technique  was  performed  by  Nisenson 
et  al  . [17].  They  confirm  the  potential  of  the  method, 
and  indicate  that  improved  phase  estimates  can  be  derived 
from  averaging  a multiplicity  of  integrated  paths,  although 
they  do  not  state  which  paths  were  used,  nor  how  many  were 
averaged.  R.V.  Stacknik  et  al  . [18],  another  group 
working  at  Itek  then  applied  the  K-T  technique  to  the 
reconstruction  of  solar  surface  features  from  speckle 
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photographs,  a test  case  providing  a very  high  SNR.  Their 
results  showed  near  d i f f r ac 1 1 on- 1 i m i ted  resolution  on  data 
obtained  from  the  24  inch  solar  vacuum  telescope  at  Kitt 
Peak,  AZ.  Both  of  these  papers  mention  the  development  of 
an  optical  processor  to  speed  the  computation  of  the 
autocorrelation  of  the  Fourier  transforms,  although  the 
processing  actually  implemented  was  entirely  digital. 

II. 3 Estimating  Phase  from  Phase  Differences  ► 

We  have  already  noted  that  the  phase  estimation  scheme  >• 

proposed  by  Knox  makes  only  limited  use  of  the  available 
phase  information,  and  that  a substantial  improvement  has 
been  reported  by  averaging  over  many  paths.  Attempting  a 
direct  extension  to  average  over  all  possible  line 

integrals  leads  to  an  unwieldy  combinatorial  problem  whose 
solution  would  be  prohibitively  expensive  for  the  large 
data  arrays  (256x256)  being  considered  here.  This  same 
problem  of  optimum  phase  estimation  from  measured  phase 
differences  also  arises  in  the  design  of  active  optical 
devices,  and  is  in  fact  a problem  in  numerical  integration. 

One  of  the  first  applications  of  1 east- squares 
estimation  to  this  phase  integration  problem  was  made  by 
Rimmer  [19]  in  evaluating  the  performance  of  a lateral 
shearing  interferometer.  Shortly  thereafter,  his  approach 
was  applied  to  the  design  of  an  active  optical  system 
intended  to  compensate  for  atmospheric  turbulence.  In  this 
system,  an  interferometer  is  used  to  measure  phase 
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differences  across  an  aperture.  These  differences  are  then 
integrated  to  form  a phase  estimate,  whose  effects  are 
removed  by  a deformable  lens  or  mirror  in  real  time.  Fried 
[20]  and  Hudgin  [21]  derived  the  same  equations  again 
almost  simultaneously.  Hudgin  assumed  that  the  s-  ution 
was  periodic,  as  when  using  the  Discrete  Fourier  transform 
(DFT),  while  Fried  did  not.  Both  researchers  derived  an 
iterative  solution  scheme,  and  demonstrated  that  for  the 
small  array  sizes  common  in  active  optics  the  rms  error  in 
the  final  phase  estimate  was  smaller  than  the  rms  error  in 
the  measured  phase  differences.  They  also  showed  that  the 
rms  phase  reconstruction  error  is  only  logarithmically 
dependent  on  the  array  size,  indicating  that  it  might  be 
applicable  to  the  large  arrays  of  speckle  imaging.  Fried 
further  reported  that  to  assure  convergence  of  his 
iterative  solutions,  he  increased  the  number  of  iterations 
almost  linearly  with  array  size.  Based  on  Hudgin's  work,  a 
group  at  Itek  headed  by  Hardy  [22]  built  an  analog  computer 
and  applied  it  to  a small  (about  5x4  array)  prototype 
active  optical  system,  with  good  results. 

In  a recent  paper  Hunt  [23]  derives  the  matrix  form  of 
these  equations  and  points  out  that  the  iterative  solution 
scheme  of  both  Fried  and  Hudgin  corresponds  to  Jacobi 
iteration,  the  numerical  iteration  scheme  with  the  slowest 
known  rate  of  convergence.  He  suggests  the  use  of  faster 
iterative  schemes,  such  as  Gauss-Seidel  or  successive 
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over- rel axati on  (SOR).  He  also  derives  an  analytical 
expression  for  the  rms  phase  reconstruction  error  based  on 
a periodic  structure  which  slightly  underbounds  the  results 
of  simulations  obtained  by  both  Fried  and  Hudgin.  Hunt's 
result  shows  the  same  nearly  logarithmic  growth  of  rms 
phase  error  with  array  size. 

Hunt's  suggestion  of  a faster  iterative  scheme  leads 
one  to  speculate  on  the  possible  use  of  a direct  algorithm 
for  solving  the  problem.  As  will  be  shown  in  chapter  III, 
the  least-squares  equations  have  a highly  regular  structure 
of  the  same  form  as  the  5-point  discrete  approximation  to 
the  Poisson  equation.  This  is  fortunate,  since  numerical 
analysts  have  in  the  last  decade  devised  some  surprisingly 
fast  algorithms  for  the  solution  of  Poisson's  equation. 
There  are  three  distinct  classes  of  algorithms,  all  of 
which  will  be  discussed  below. 

The  first  is  the  method  of  Cyclic  Reduction  (CR) 
devised  by  Buneman  [24].  He  noted  that  the  even-numbered 
equations  could  be  combined  to  form  an  equation  for  every 
fourth  line,  which  could  in  turn  be  reduced  to  an  equation 
for  every  eighth  line,  etc.  There  is  finally  only  one 
equation  remaining,  which  is  easily  solved.  The  solution 
on  the  odd-numbered  lines  is  obtained  from  the  reduced 
equations  corresponding  to  the  successive  levels  of 
reduction.  The  method  is  particularly  easy  to  code,  and  is 
numerically  stable.  It  is  not  the  fastest  method,  however. 
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The  third  type  of  algorithm  is  that  of  Giwa  [25],  and 
is  an  adaptation  of  an  analytical  va r i a t i on- of- parameter 
method  for  solving  second-order  partial  differential 
equations.  It  is  very  general,  highly  stable,  and  does  not 
restrict  the  dimensionality  of  the  problem,  as  do  the 
Fourier  methods.  It  is  less  memory- ef f i c i ent  and  slower 
than  the  other  methods,  however. 

The  fastest  known  algorithm  is  that  of  Hockney  r?6l, 
which  combines  both  FA  and  CR.  It  is  the  most  difficult  to 
code,  and  must  be  carefully  used  to  avoid  overflow 
problems.  In  his  excellent  review  paper,  Hockney  [271 
compares  his  own  algorithm  with  CR  and  the  best  iterative 
schemes,  including  SOR.  He  indicates  that  his  FACR  is  the 
fastest,  that  FA  is  about  507  slower,  and  CR  is  1001- 
slower.  He  also  shows  that  to  guarantee  a maximum  error  of 
one  part  in  10F6,  about  4.1?  iterations  aVe  required  of  SOR 
on  an  array  of  128x128.  Since  the  direct  methods  achieve 
the  same  accuracy  in  a time  equal  to  that  required  by  about 
six  SOR  iterations,  there  seems  to  be  no  reason  to  use  an 
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iterative  method . 

II. 4 Implicit  Phase  Estimation 

Implicit  phase  estimation  is  the  name  I have  given  to 
the  technique  of  estimating  the  phase  spectrum  based  on  a 
priori  knowledge  of  the  signal  properties  and  the  magnitude 
spectrum  alone,  rather  than  on  direct  measurements  of  the 
phase  spectrum  or  its  derivative.  The  approach  has  appeal 
because  it  is  usually  much  easier  to  measure  accurately 
magnitude  spectra,  as  in  the  case  of  electron  microscopy, 
x-ray  crystallography,  and  image  restoration  [30].  It  has 
been  explored  by  many  authors  over  the  past  two  decades 
[31-36],  and  remains  an  area  of  active  interest. 

For  one- dimens i onal  signals,  the  problem  may  be  easily 
understood  in  the  framework  of  linear  system  theory.  If 
the  duration  of  such  a signal  is  finite,  then  its  Fourier 
transform  is  an  analytic  function.  Thus,  the  real  and 
imaginary  parts  of  the  transform  are  related  by  the 
Cauc hy -R i emann  equations.  Other  analytic  spectra  may  be 
formed  however  by  reflecting  any  transform  zero  to  its 
conjugate  reciprocal  location.  This  will  modify  the  phase 
spectrum  while  leaving  the  magnitude  spectrum  unchanged. 
Wolf  [33]  derived  the  allowable  forms  of  the  resulting 
phase  spectra  assuming  the  object  to  be  real,  thus 
demonstrating  that  the  object  reconstruction  is  not  unique. 
However,  in  two  dimensions  the  situation  is  not  as  easy  to 
describe.  This  is  because  analysis  in  terms  of  transform 
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poles  and  zeroes  is  not  generally  possible.  The  lack  of  a 
fundamental  theorem  of  algebra  for  functions  of  two 
independent  variables  makes  spectral  fac tor i za ti on 
techniques  inapplicable. 

Nevertheless  a number  of  schemes  for  implicit  phase 

estimation  have  been  proposed.  Kohler  and  Mandel  T32] 

compare  two  methods  on  one  dimensional  spectral  data.  The 

first  involves  the  use  of  a reference  point  source,  the 

second  a numerical  integration  of  a Hilbert  transform  f 

relationship  when  the  zero  locations  in  the  z plane  are  >• 

l 

approximately  known.  They  indicate  that  the  second  method 
requires  very  accurate  knowledge  of  the  magnitude  spectrum, 
and  their  method  does  not  ?eem  to  have  been  applied  to  real 
astronomical  data.  p 

l s 

In  electron  microscopy  a slightly  different  situation 

' r 

obtains  in  that  the  magnitude  spectrum  is  known  in  both  the 
image  and  diffraction  planes  of  an  imaging  system  while  the 
phase  is  not  known  in  either  plane.  It  is  well  known  that 

1 J 

the  signals  in  these  two  planes  are  Fourier  transform 

j 

pairs.  Gerchberg  and  Saxton  [37]  propose  an  iterative 
algorithm  for  determining  the  phase  in  both  planes.  The 
algorithm  assumes  an  initial  random  phase  distribution  in 
one  of  the  domains.  The  Fourier  transform  is  taken,  and 
the  resulting  magnitude  spectrum  is  corrected  to  agree  with 
the  measured  magnitude  spectrum  in  that  domain.  The 
inverse  transform  is  then  taken  and  its  magnitude  spectrum 
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is  similarly  corrected.  The  entire  process  iterates  until 
convergence,  which  may  be  quite  slow.  The  uniqueness  of 
the  resulting  solution  has  been  examined  in  [38],  where  the 
authors  conclude  that  the  reconstructed  phase  is  unique 
provided  that  the  magnitude  spectrum  does  not  have  even 
symmetry  and  is  non-zero,  and  that  the  resulting  frequency 
spectrum  is  analytic.  If  the  magnitude  spectrum  is  even, 
there  may  exist  at  most  one  other  analytic  function 
satisfying  the  constraints.  Methods  to  speed  the  algorithm 
convergence  have  been  described  by  both  Gassman  [39]  and 
Feinup  [40].  The  latter  has  also  experimented  with  this 
algorithm  when  only  the  magnitude  spectrum  is  known  in  the 
frequency  domain,  and  only  the  region  of  support  is  known 
in  the  signal  domain.  While  solution  non-uniqueness  is  a 
serious  problem  for  one  dimensional  signals,  he  obtains 
unique  solutions  on  simulations  using  complicated  two 
dimensional  images.  It  is  not  presently  known  why  the 
two-dimensional  results  seem  to  be  unique.  It  is  known 
[41]  that  an  initial  (even  highly  inaccurate)  phase 
estimate  can  often  be  used  to  distinguish  the  true 
solution,  as  well  as  to  greatly  speed  convergence  of  the 
algorithm.  Feinup  proposes  the  use  of  this  method  to  solve 
the  phase  problem  of  speckle  interferometry. 

II. 5 Super- resol ut i on 
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The  extrapolation  of  the  frequency  spectrum  of  a 
lowpass  filtered  signal  is  called  super- resol ut ion , because 


its  effect  is  to  increase  the  apparent  detail  or  resolution 
of  the  reconstructed  signal.  A number  of  techniques  have 
been  proposed  to  extrapolate  the  low-frequency  spectrum. 
Harris  [42]  proposed  approximating  the  true  spectrum  by 
forming  a weighted  sum  of  sine  functions,  where  the  weights 
are  the  coefficients  of  the  discrete  Fourier  series  of  the 
signal.  The  extended  spectrum  appears  to  be  sensitive  to 
noise  in  the  data,  however.  Barnes  [43]  expands  the  given 
portion  of  the  spectrum  in  terms  of  prolate  spheroidal  wave 
functions  [44].  It  is  shown  by  Rushforth  and  Harris  [45] 
that  this  method  is  also  highly  sensitive  to  measurement 
error,  and  that  achieving  resolution  much  beyond  the 
diffraction  limit  with  data  of  realistic  accuracy  Is 
doubtful.  Frieden  [46]  describes  the  'maximum  entropy' 
approach,  which  is  probabilistic  in  nature.  Although 
substantial  computation  is  required,  good  results  have  been 
obtained.  Gerchberg  [47]  describes  an  iterative  scheme 
which  he  calls  error-energy  reduction.  He  shows  that  the 
algorithm  is  linear  and  that  given  noise-free  data  the 
solution  converges  to  the  true  spectrum.  Convergence  can 
be  quite  slow.  Gerchberg  notes  that  error  in  the 
extrapolated  portion  of  the  spectrum  due  to  measurement 
noise  is  increased  at  every  iteration,  thereby  limiting  the 
accuracy  of  reconstruction  possible  on  real  data.  An 
analysis  of  the  algorithm  by  Youl a [48]  casts  the  problem 
in  the  form  of  alternating  orthogonal  projections  on  linear 
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subspaces  of  a parent  Hilbert  space.  His  analysis  also 
shows  the  ill-conditioned  nature  of  the  algorithm.  In 
particular  he  notes  that  there  is  an  optimal  number  of 
iterations  after  which  the  error  in  the  continued  spectrum 
increases  rather  than  decreases,  and  that  this  number 
cannot  in  general  be  known  a priori.  Still,  this  technique 
seems  to  be  the  most  successful  to  date,  and  an  essentially 
identical  algorithm  has  been  used  by  Papoulis  [49]  to 
extrapolate  a truncated  band-limited  time  waveform. 
However,  a new  algorithm  based  on  Youla’s  work  has  recently 
been  described  by  Cadzow  [50],  and  it  may  prove  superior 
although  it  has  not  yet  been  tested  on  noisy  data. 


CHAPTER  III 

LEAST  SQUARES  PHASE  ESTIMATION 

1 1 1 . 1 Introduction 

In  this  chapter  I describe  a method  of  estimating  the 
phase  spectrum  of  an  object  whose  image  is  both  blurred  and 
noisy.  The  method  is  very  similar  to  power  spectral 
density  (PSD)  estimation,  and  should  have  application  in 
the  same  situations.  The  basis  of  the  method  is  the  use  of 
the  Knox-Thompson  estimate  of  phase  differences.  Section 

111. 2 contains  a brief  explanation  of  the  fundamentals  of 
speckle  interferometry  and  the  Knox-Thompson  estimator. 
However,  for  detailed  analysis,  the  work  of  Knox  [15] 
should  be  consulted.  In  section  III. 3 the  least-square 
phase  estimate  is  derived,  as  well  as  a fast  direct  Fast 
Fourier  Transform  (FFT)  based  algorithm  to  solve  the 
resulting  equations.  Section  III. 4 will  discuss  the 
algorithm  performance  in  terms  of  speed  and  integration 
error,  and  section  III. 5 will  compare  the  least-square 
phase  estimate  with  a two-path  line  integral  phase  estimate 
for  various  SNR  on  simulated  data.  The  effects  of  spatial 
windows  on  the  estimates  will  also  be  noted,  as  will  the 
use  of  phase  difference  information  derived  from  the  second 
shift  term  of  the  transform  autocorrelation.  Further 
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refinement  of  the  phase  estimates  with  a Gere hberg - Sax  ton 
algorithm  is  demonstrated  in  section  f> . This  section  also 
contains  a proof  of  the  convergence  of  the  algorithm  in 
this  application,  and  a proof  of  the  possible 
non- conv ergence  of  a similar  refinement  operation  on  the 
magnitude  spectrum. 

III. 2 Review  of  Speckle  Imaging 

As  mentioned  in  the  introduction,  atmospheric 
turbulence  causes  a variation  in  the  refractive  index  of 
the  optical  path.  This  is  modeled  by  a point  spread 
function  (PSF)  which  varies  in  both  the  time  and  spatial 
coordinates  of  the  imaging  system.  The  resulting  blurring 
process  must  be  modeled  as  a complicated  superposition 
integral.  However,  if  the  imaging  is  restricted  to 
sufficiently  small  solid  angles,  the  PSF  is  nearly 
spatially  shi f t- 1 nv ari ant . Furthermore,  for  short-exposure 
images,  Labeyrie  showed  that  the  physical  process  is 
approximately  convolutional.  The  exposure  times,  typically 
5 to  10  milliseconds,  are  so  short  that  the  atmosphere  can 
be  regarded  as  a fixed  optical  medium,  whose  index  of 
refraction  remains  spatially  random.  The  primary  effect  of 
the  atmosphere  on  each  short- ex po sure  photograph  is  to 
Introduce  a random  phase  component  to  the  spectrum,  causing 
constructive  and  destructive  interference  of  the  light 
waves  to  occur  randomly  in  the  image  plane.  As  a result, 
the  photographs  are  composed  of  many  small  patches  of 
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light,  or  speckles,  hence  the  name  speckle  imaging. 

Consider  an  incoherent  imaging  system  with  intensity 
distribution  i(x),  where  x is  a 2-dimensional  position 
vector  describing  a location  in  the  focal  plane  of  the 
system.  Over  a shi f t- i nv ar i ant  ( isopl anatic)  region  i ( x ) 
is  related  to  the  object  intensity  distribution  i Q( x ) by 
convolution,  as 

i ( X ) = J i o ( a ) s ( X+  f a / z ) da  , ( 1 ) 
where  s ( a ) represents  the  instantaneous  point- spread 
function  (PSF)  of  the  atmosphere- tel escope  combination,  f 
is  the  effective  focal  length  of  the  telescope,  and  z is 
the  distance  between  the  object  and  telescope.  The  Fourier 
transformation  of  this  equation  is 

I(u)«I*(fu/z)S(u)  , (2) 
where  capital  letters  denote  the  Fourier  transform  of  the 
lower-case  quantity.  Because  convolution  is  a linear 
operation,  a long-exposure  image  can  be  equivalently 
represented  by  the  average  of  the  Fourier  transforms  of 
many  short- expo  sure  images  in  the  Fourier  domain,  as 

< I ( u)>= I *( fu/z) <S ( u) > , (3) 
where  the  angle  brackets  denote  averaging.  It  has  been 
shown  both  theoretically  and  experimentally  that  <I(u)>  is 
a narrow  function,  i.e.  has  a low  cutoff  frequency.  The 
heart  of  Labeyrie's  technique  is  the  fact  that  under  the 
conditions  required  to  produce  speckle  images  the 
mean-square  OTF  ( < I S ( u ) I >)  has  non-zero  magnitude  out  to 
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the  di f f r ac t i on- 1 i m i t of  the  telescope.  The  phase  of  the 
image  transforms  has  suffered  a random  phase  perturbation, 
however.  If  one  attempts  to  average  directly  many 
short- ex posure  images  or  their  transforms  the  random  phase, 
which  is  most  pronounced  at  the  higher  frequencies,  causes 
a cancellation  of  the  high  spatial  frequency  energy.  This 
cancellation  can  be  avoided  if  one  averages  instead  the 
squared  modulus  of  the  transforms,  as  first  shown  in  this 
context  by  Labeyrie  [ 1 J . This  averaging  essentially 
ignores  the  phase  contribution,  whatever  it  may  be.  Thus, 
one  forms 

< I ( u ) I ( u ) > = I Q ( f u / z ) IQ( fu/z)<S(u)S*(u)>  (4) 

or  equivalently 

<|l(u)|2>=|lo(fu/z)|2<|S(u)|2>.  (5) 


The  squared  modulus  of  the  OTF  is  then  approximated  by 
observations  of  an  unresolvable  point  star,  and  the  PSD  of 
the  object  is  calculated  by  simple  inverse  filtering,  as 


l!o 

( fu/z)  I 

2=<|  I ( 

u) | 2>/<| S(u)  | 2>  . 

(6) 

It  is 

wel  1 

known 

that 

inverse  f i 1 ter i ng 

i s 

often 

inadequate 

since 

it  tends 

to  amplify  noise 

a t 

those 

frequenc i es 

where 

there 

is  little  or  no  energy  in 

the 

OTF  . 

It  has  not  been  found  necessary  in  our  work  to  use  other 
filtering  strategies,  however,  for  the  following  reasons. 
First,  because  it  has  essentially  a triangular  or  Bartlett 
frequency  response,  the  PSD  of  the  OTF  has  appreciable 
energy  in  all  spatial  frequencies  out  to  the  absolute 
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frequency  cutoff  imposed  by  the  d i f f rac ti on- 1 imi t of  the 
telescope.  Second,  there  is  non- neg 1 ig ib 1 e noise  energy  in 
our  measurements,  and  averaging  together  many  noisy 
observations  tends  to  put  some  noise  energy  in  all  the 
frequencies.  If  enough  images  are  averaged,  the  inverse 
filtering  operation  is  stable. 

PSD  estimation  is  presently  a well-established 
technique  in  many  fields,  and  is  described  in  a number  of 
texts  [56-58].  By  the  We i ner- Kh i n tc i ne  theorem  the  PSD  and 
the  au tocor rel a t i on  of  a given  function  are  a Fourier 
transform  pair.  Therefore  the  Labeyrie  estimate  is 
equivalently  represented  as  the  object  autocorrelation. 
This  is  sufficient  for  c ha rac ter i z i ng  centrosymmetric 
objects  and  for  measuring  interstar  distances  for  binary 
stars.  However,  without  knowledge  of  the  object  phase,  no 
object  reconstruction  can  be  made  in  the  general  case. 

The  PSD  can  also  be  described  as  the  zero-shift  term 
in  the  autocorrelation  of  the  Fourier  transform.  The 


Knox- Thompson  estimator  is  similar  to  that  of  Labeyrie  in 
that  It  uses  terms  (other  than  the  zero-shift  term)  from 
the  statistical  autocorrelation  of  the  image  transform, 
this  time  to  estimate  phase  differences.  Form  the  product 
< I ( u 1 ) I ( U2) > = I 0 ( fu 1 / z ) I 0 ( f u2/ z) < S ( u1 ) S (u2)>*  (7) 

If  I o ( u ) is  represented  in  polar  form 

IQ( fu/z)  = f IQ( fu/z) | exp(j©( fu/z) ) (8) 

then  it  is  easy  to  see  that  the  phase  differences  can  be 
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recovered  from  the  expression 
<1  (u,  )I  (u-)> | <S(u,  )S*(u„)> | 

r < I ( u , ) 1 « ( u2 ) > |~s 1 5, ) S*  < )T  ■e«PfJe(fu2/z)-j.(f„]/2)).  (9) 
The  limitations  of  this  approach  are  derived  theoretically 
by  Knox  in  his  thesis,  assuming  that  the  effects  of 
atmospheric  turbulence  are  adequately  approximated  by  a 
random  phase  model  with  large  variance,  i.e.  greater  than 
one  wavelength.  (While  it  is  known  that  some  random 
amplitude  variations  are  also  introduced  by  the  atmosphere, 
it  is  generally  agreed  that  the  phase  variations  have  a 
much  greater  effect  on  image  formation).  Large-scale  phase 
variations  are  removed  by  recentering  the  speckle  images 
about  their  centroids,  so  that  the  remaining  phase 
variation  is  V ( x ) = ex  p { j $(  x ) } , where  <t>(x)  is  a stationary 
Gaussian  random  process.  Knox  shows  that  the  mean  OTF  can 
be  written  as  a product  of  functions  describing  separately 
the  effects  of  the  atmosphere  and  the  telescope,  and  that 
the  atmospheric  'seeing'  defines  the  cutoff  frequency  of 
the  mean  OTF.  His  asymptotic  evaluation  of  the 
autocorrelation  of  the  OTF  indicates  the  presence  of  a bias 
term  with  non-zero  value  out  to  the  d i f f r ac t i on- 1 i m i t if 
both  u1  and  u^  are  zero,  or  if  they  are  both  large  and 
their  difference  is  small.  The  presence  of  this  bias  term 
for  u i and  u2  small  is  assumed  on  physical  grounds,  since 
it  is  present  in  the  limit  whether  u^  and  u^  grow  large  or 
approach  zero.  (This  intuition  seems  to  be  justified  based 
on  the  results  of  simulations).  Knox  also  shows  that  the 
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technique  is  quite  insensitive  to  telescope  aberrations. 
He  notes  that  it  is  sensitive  to  noise  and  that  for  best 
results  the  distance  between  u -j  and  U2  should  be  kept 
small,  so  as  to  maximize  the  SNR. 


III. 3 Derivation  of  the  Lea s t- Squa res  Phase  Integrator 

Since  phase  estimation  appears  to  be  the  primary 

I 

difficulty  in  imaging  through  the  atmosphere,  one  would 
like  it  to  be  robust,  make  maximum  use  of  the  data  (in  some 
sense),  and  be  cheap  computationally.  Phase  estimation  by 
averaging  a few  line  integrals  is  easy  to  compute,  but 
obviously  makes  only  limited  use  of  the  phase  difference 
data,  and  as  a result  produces  phase  estimates  of  greater 
variance  than  might  be  desired.  It  also  requires  the  use 
of  arithmetic  performed  modulo  2it  (in  order  that  the 
averaging  of  different  paths  be  consistent),  and  delivers 
the  principal  value  of  the  phase,  rather  than  the  phase 
itself.  A 1 east- squares  approach  has  the  well-known 
advantage  of  satisfying  a tractable  objective  function.  In 
addition,  it  does  use  the  data  in  a more  'general'  way,  and 
in  practice  leads  to  a robust  estimator  in  noise  whose 
statistics  approximate  (sometimes  quite  loosely)  the 
Gaussian.  It  requires  only  conventional  arithmetic  (with 
the  proviso  that  the  Fourier  tr  nsform  be  sampled  closely 
enough  so  that  the  phase  differences  are  themselves  less 
than  2"),  and  it  returns  an  'unwrapped',  or  non- pr i nc i pa  1 
value  phase  estimate.  The  difficulty  to  date  has  been  the 
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huge  dimensionality  of  the  problem.  As  in  many  other  image 
processing  problems,  an  image  represented  by  an  array  of 
256  by  256  discrete  points  expands  to  a linear  system  of 
equations  having  65536  equations  in  as  many  unknown s . 
♦ Iterative  numerical  schemes  have  previously  been  proposed 

to  solve  the  problem  [20-23].  However,  equally  accurate 
and  much  faster  direct  algorithms  are  available,  as  noted 
» i n Chapter  1 1 . 3 . 

We  have  chosen  an  FA  algorithm  for  the  following 
reasons.  First,  since  we  are  restricted  to  digital 
« processing  of  the  speckle  images,  our  use  of  the  FFT 

already  requires  the  appropriate  dimensioning  of  image 
arrays,  and  no  further  restriction  is  encountered.  Second, 
the  algorithm  is  nearly  as  fast  as  FACR,  but  is  much  less 
difficult  to  code.  Third,  it  has  a neat  modular  structure 
that  allows  its  use  on  minicomputers  with  very  small 
core-memory  stores.  In  the  following  we  will  present  the 
derivation  of  the  least-squares  equations  and  two  solution 
algorithms.  The  first  assumes  non-periodic  boundary 
conditions,  and  so  is  an  approximation.  It  might  prove 
useful  should  optical  Fourier  transform  techniques  be  used 
in  data  processing.  For  our  work  the  periodicity  of  the 
discrete  Fourier  transform  (DFT)  implies  periodic  boundary 
conditions,  however,  and  in  this  rase  the  second  algorithm 
is  exact. 
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III. 3.1  The  Normal  equations 

Consider  a rectangular  sampled  image  transform  of 
dimension  M*N,  and  denote  the  phase  associated  with  point 
(i,j)  by  $.jj.  We  define  noisy  phase  differences  in  the  two 
directions  by 

r,  . . i . + e . . 

1J  1J  1-1 ,J  13 

and 


C.  .=«!>.  .—<#>.  . ,+r).. 

1 J 1 J 1 ,3-1  1 J 

where  and  n.^  represent  deviations  from  the  true  phase 

differences  associated  with  the  object.  Our  objective  is 

A 

to  obtain  estimates  ^ j at  each  point  in  the  image 
transform  by  integrating  the  noisy  phase  differences.  The 
problem  is  o ver- de term i ned , and  we  use  a 1 ea s t- squa re s 
approach  to  phase  estimation. 

We  assume  that  the  noise  statistics  are 
shi f t- i nv ar i ant  and  define  the  sum  of  squared  errors  by 

. 2 


M N 

-.1. 
i=l  J=1 


T . .+$ . , .J^+fc.  . . )2]  (10) 

13  i J 1-1 .3  13  13  i ,3-1 


initially  assuming  that  is  zero  outside  the  image 

array.  Differentiating  (10)  with  respect  to  4>  • • and 

* \J 

equating  the  result  to  zero  leads  to  the  following 
equations: 

4(i>.  .-j>.  . , .,=r,,-r..,  .+c..-c.  ... 

13  1-1.3  1 + 1, j 1,3-1  i.J  + 1 13  i + l, J 1J  1,3  + 1 

2< i <M- 1 , 2<j<N-l  (11) 

3*ij-*2j_<h  ,j-r*i  j+rci  j"ci  ,j+rr2,j 

i = l , 2< j <N-  1 (12) 
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Equation  (11)  is  valid  at  an  interior  point  of  the  image 
transform  array,  (12)-(15)  at  the  edges,  and  (16)-(19)  at 
the  corners.  However,  if  the  phase  array  is  assumed  to  be 
periodic,  as  it  is  in  our  work  with  the  speckle  images 
because  of  the  periodicity  of  the  DFT,  then  (11)  is  valid 
for  all  points  l<i<M,  l<j<n,  provided  i-1  and  j-1  are  taken 
modulo  M and  N. 

The  equations  and  the  algorithms  for  solving  them  in 
the  periodic  and  aperiodic  cases  differ  slightly,  and  the 
algorithms  will  be  described  separately.  In  both  cases, 
the  equations  to  be  solved  are  analogous  to  those  which 
arise  in  solving  Poisson's  equation  on  a rectangle  using 
the  standard  five-point  finite  difference  approximation 
[24-28].  In  particular,  we  discuss  two  variants  of  an 
FFT-based  algorithm  described  in  [28]. 


III. 3. 2 The  Aperiodic  case 
We  first  define 

virrij-ri*i,o+circi,o*i  (zo) 

where  r..j  and  c^j  are  noisy  phase  differences.  We  then 
consider  the  v^j  to  be  arranged  in  a block  column  vector, 
with  v..  being  the  jth  element  of  the  ith  block.  The  phase 

* J 

elements  are  ordered  in  an  identical  manner.  For  large 

* J 

M and  M (the  usual  case  in  image  processing),  (11)-(19)  can 
be  approximated  by  the  system  of  matrix  equations 

A*  = v (21) 


s: 

4 


l J 


' 

I 1 


where  the  coefficient  matrix  A has  the  block  form 
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wi  th 


(22) 


(23) 


A Is  block  tri-diagonal  of  block  dimension  M*N,  I is  the 
N*N  identity  matrix,  and  A0  is  an  N*N  tri-diagonal  matrix. 
The  block  matrix  can  be  written  in  the  alternative 

form 


Vr'r'i 

(24) 

and 

(25) 

(26) 

Since  Aq  is  symmetric,  there  exists  an  orthogonal  matrix  Q 
with  Q ' Q = I , where  the  prime  denotes  matrix  transpose,  and 


Q'AoQ=D0.  (27) 

As  is  well  known,  the  matrix  Q has  as  its  columns  the 
normalized  eigenvectors  of  AQ , so 

q1js  (|f^T^Z  sin^fUT^  1»j  = l»2,...N  (28) 
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Mote  that  for  this  case  0=0 ' • D0  is  a diagonal  matrix 
whose  diagonal  elements  are  the  eigenvalues  of  A0 : 

Doit=4-2  c°s<NTr>  <29> 

Using  (27),  we  can  write  (24)-(26)  in  the  form 


where 


A1  ternatively , 


oVVvi 

(30) 

♦Vn-n+r*! 

(31  ) 

M-1+0o4M=vM 

(32) 

(33) 

V«,vi 

(34) 

Do  -1 


-I  Dq  -I 


Permuting  the  rows  and  columns  of  this  coefficient  matrix 
by  grouping  the  1th  equation  in  each  block  together  results 
in  the  alternative  coefficient  matrix 

r«,  /~\  l 
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In  other  words,  the  original  block  tri-diagonal  system  can 
be  decomposed  into  M decoupled  tri-diagonal  systems,  each 
of  which  can  be  quickly  solved  using  Gauss  elimination. 
Having  solved  these  decoupled  systems,  we  then  undo  the 
permutation  described  above  and  premultiply  the  result  by  Q 
to  obtain  the  solution  to  the  original  system. 

We  see  from  (28)  the  the  diagonal ization  process 
essentially  amounts  to  calculating  a DFT.  If  N + l and  M+l 
are  highly  composite,  in  particular  if  they  are  powers  of 
2,  this  computation  can  be  performed  rapidly  using  the  fast 
Fourier  transform  (FFT)  algorithm.  To  accomplish  this,  we 
first  define  for  each  k 


Re(tt)= 


kt 


e = 0 

t = l ,2,...N 
< =N  + 1 ....  2 (N  + l )-l 


and 


(38) 


I m f t { > = 0 


e = 0,l ,...2(N  + 1 )-l 


(39) 
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Then 

o 1/2 

vke=(RTT)  Im{DFT(t*)}  *■  = ! ,2 , . . .N  (40) 


where  the  DFT  is  defined  by 


t 

DFT ( 1 1 ) 

2 (N  + l )- 

■ l 

1 = 0 

teexp[2TTnry] 

(41  ) 

Storage  requirements 

and  computation  time  can 

be 

0 

further 

reduced  by 

taking 

advantage  of  the  even  and 

odd 

symmetry 

properties 

of  the 

DFT.  Because  they  are 

wel  1 

known,  we  do  not  discuss  them  here,  although  our  programs 
♦ do  incorporate  them. 


III. 3. 3 The  Periodic  case 

The  strategy  in  the  periodic  case  is  similar  to  that 
outlined  above.  In  this  case,  the  coefficient  matrix  is 
the  block  circulant  matrix 


C -I  0 ...  -I 
o 

-I  C -I  ...  0 

o 

• 

-1  0...-I  C0 


(42) 


where 


-1  0 ...  -1 

4-10  ...  0 


0 ...  -1  4 


(43) 


1 


i* 


A 


A 
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For  this  case,  Q = Q ‘ . Q has  as  two  of  its  columns 

( 1 / N )V  2 ( 1 , 1 1)'  and  (1/N)V2(1,-1,1,...,-1)'  , with 

corresponding  simple  eigenvalues  2 and  6,  respectively.  In 
addition,  there  are  (N-2/2)  double  eigenvalues 

ya  = 4-2  cos(2ttS,/N)  4 = 1 , . . . (N-2)/2  (44) 

each  associated  with  a pair  of  eigenvectors 

"(n)  s i n ( — n — ) (45) 


n(l)_,2)]/2  r ,k42tr, 
nt/  (n'  COS(  jj  ) 


Di ag onal i zation  is  again  carried  out  using  the  FFT. 
In  this  case,  we  require  N to  be  highly  composite, 
preferably  a power  of  two.  After  di agonal ization  and 
permutation  as  in  the  aperiodic  algorithm,  the  decoupled 
block  coefficient  matrices  become 


-1 


-1 


0 . . . 0 


V 0 -1 


-1  0 . . . 0 


where 

Y^  = 4-2  cos(2tt£/N),  4 = 0,1,,,,  N/2  (48) 
As  described  above  for  the  aperiodic  case,  these  decoupled 
systems  can  now  be  solved  separately.  The  permutation  and 


d i agona 1 i za t i on  operations  are  then  reversed  to  obtain  the 
solution  to  the  original  system  of  equations. 

There  are  some  important  differences  between  the 
periodic  and  aperiodic  cases.  First,  rp  is  singular, 
reflecting  the  singularity  of  C.  If  the  first  row  and 
column  of  rQ  are  deleted,  a solution  is  obtained  which 
satisfies  the  original  equations  and  which  differs  from  any 
other  solution  by  an  additive  constant.  This  ambiguity  can 
be  resolved  in  our  case  by  noting  that  the  phase  at  the 
origin  must  be  zero. 

Second,  the  remaining  blocks  of  coefficients  are  not 
tri-diagonal  since  they  have  non-zero  corner  elements. 
This  situation  can  be  dealt  with  using  the  Sherman-Morr i son 
formula  [28],  which  states  that  if  A=B  + uv',  where  u and  v 
are  column  vectors,  then 

A” 1 = B" 1 +B_1U(  I + v'B'1  u ) " 1 v 1 B~ 1 
If  v contains  only  one  nonzero  element,  the  correction  is 
quickly  obtained.  Since  most  of  the  eigenvalues  are  of 
order  two,  two  blocks  can  usually  be  solved  at  once, 
thereby  compensating  somewhat  for  the  additional 

computation  involved  in  the  Sherman-Morr i son  formula.  A 
further  saving  is  realized  from  the  increased  efficiency  of 
the  FFT  in  diagonalizing  the  periodic  structure. 

III. 4 Algorithm  performance 

A number  of  tests  were  run  with  both  the  periodic  and 
aperiodic  integration  algorithms  to  determine  their 
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accuracy.  The  first  tests  related  to  the  propagation  of 
numerical  errors  irrespective  of  any  measurement  errors  in 
the  data.  The  other  tests  examined  the  propagation  of 
measurement  error  in  the  integrated  solution. 

The  first  tests  were  as  follows.  An  a ppropr i a tel y 
dimensioned  file  was  filled  with  Gaussian  distributed 
random  noise  of  standard  deviation  1.  The  discrete  Poisson 
differential  operator  produced  another  file,  which  upon 
integration  was  compared  to  the  first.  Typical  results  for 
the  algorithms  are  summarized  in  Tables  1 and  2. 

Integration  time  for  the  255*255  aperiodic  array  was 
about  200  seconds  on  our  single  user  PDP-10  computer.  (It 
has  floating  point  add  and  multiply  times  of  about  5 us  and 
11  us,  respectively).  Approximately  175  seconds  were  spent 
on  the  FFT  compu ta t i ons . Because  the  periodic  algorithm  is 
better  suited  to  an  FFT  Implementation,  the  256*256  problem 
required  only  135  seconds,  of  which  approximately  100 
seconds  were  spent  on  the  FFT.  Corrections  by  the 
Sherman-Morr i son  formula  required  only  about  10  seconds. 

Previous  researchers  have  indicated  the  stability  of 
the  growth  of  the  mean  squared  error  (MSE)  in  the 
Integrated  solution  due  to  errors  in  estimating  the  phase 
differences  when  the  1 east- squares  approach  is  used.  Both 
Fried  and  Hudgin  estimated  the  growth  as  logarithmic  with 
Increasing  array  size.  Based  on  simulations  with  small 
arrays.  Fried  found  that  a good  fit  to  the  mean  square 


Table  1 


Error  in  the  estimated  phase  due  to  numerical 
errors  when  using  the  aperiodic  integrated 
algorithm. 
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Table  2 
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Array  Size 

Maximum  error 

SSE  Norm 

15x15 

. 33x 1 O-6 

- 1 3 

.14x10  IJ 

31x31 

. 48x1 O-6 

- 1 3 

.18x10  1 

63x63 

. 9x 1 0"6 

- 1 2 

.18x10 

127x127 

. 8x1 0-5 

- 1 2 

.37x10 

255x255 

. 23x 1 0~  5 

. 80x 1 0” 1 1 

Error  in  the  estimated  phase  due  to  numerical 
errors  when  using  the  periodic  integrated 


algori thm. 

Array  Size 

Maximum  Error 

SSE  Norm 

16x16 

. 2 1 x 1 0”6 

.73xl0"14 

32x32 

.66xl0"6 

.98x1  O'1 3 

64x64 

. 98x 1 O-6 

- 1 2 

.83x10  l£ 

128x128 

.20xl0'5 

- 1 2 

.41x10 

256x256 

. 48x 1 O'  5 

.46x1  O'11 

J 
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reconstruction  error  could  be  made  with  the  function 

e =Opd[.  3205  1 n ( N ) + .6558]  (49) 

2 

where  opd  is  the  variance  of  the  errors  made  in  estimating 
the  phase  differences.  Similarly,  Hudgin  obtained  the 
function 

e =o2  [.103  1 n ( N ) + .561]  (50) 

pd 

Hunt  derives  an  expression  for  the  error  based  on  periodic 
boundary  conditions,  which  he  gives  as 

e =^apd^2^  I [2  si  n(  i n/N  ) /4-2cos(  2 iri /N  ) -2cos(  2 /N  ))  ] 

+ [2si  n(  j n/N  )/ (l-2cos(  2 iri /N  )-2cos(  2 irj/N  ) ) ]2} . (51) 
He  shows  that  the  mean  square  error  in  the  integrated  phase 
estimate  predicted  by  his  trigonometric  expression  has  the 
form  of  the  empirically  derived  logarithmic  expressions  of 
Fried  and  Hudgin.  He  also  notes  that  upon  plotting  the 
results  there  is  a distinct  offset  between  his  predictions 
and  the  results  of  the  other  two.  This  offset  is  ascribed 
to  the  different  ways  in  which  their  algorithms  actually 
handled  the  boundary  conditions.  A plot  of  these  three 
curves  together  with  typical  solutions  from  our  periodic 
algorithm  is  found  in  Fig.  1. 

The  stability  of  the  reconstructions  with  increasing 
array  size  is  encouraging.  Note  also  that  the  MSE  of  the 
reconstructed  phase  estimates  is  actually  less  than  the 
variance  of  the  phase  difference  measurement  error  for 
array  sizes  up  to  256*256.  Because  of  the  periodicity  of 
the  DFT  Hunt's  assumption  of  periodic  boundary  conditions 
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Fig.  1. 


Growth  of  the  ratio  of  mean  square  error  in  the 
integrated  phase  estimate  to  the  variance  of  the  error  in 
measuring  phase  differences  as  a function  of  array  size  N. 


is  exactly  satisfied,  and  our  results  agree  quite  well  with 
his  analytically  derived  result. 

III. 5 Phase  Estimation  on  Simulated  Data 

Simulations  were  used  extensively  to  evaluate  the 
effectiveness  of  the  various  estimation  procedures 
discussed  in  this  research.  The  method  of  generating 
synthetic  data  is  described  below.  Note  that  the  resulting 
synthetic  data  base  is  common  to  all  the  simulations  in 
this  dissertation. 

A synthetic  binary  star  was  constructed  of  a separable 
triangle  of  height  l and  a shifted  and  scaled  version  of 


Itself;  it 

is  shown 

i n 

Fig. 

2a  . It  i s def i ned 

on  a 

6 4*64 

array  and 

has  a 

reg 1 on 

of  support  of  18 

pixel s . 

Its 

assocl ated 

magnitude 

and 

phase  spectra  are 

shown 

i n 

Figs.  2b 

and  2c, 

respectively.  (Note 

that 

the 

illustrations  of  the  frequency  spectrum  are  represented  in 
the  base  band  with  the  zero  frequency  in  the  upper  left 
hand  corner.  They  are  also  transposed  relative  to  the 
intensity  object,  to  save  processing  time.  This  convention 
will  hold  throughout  the  dissertation).  This  object  is 
blurred  by  a digital  process  Imitating  the  effects  of 
atmospheric  turbulence.  We  assume  that  the  distortions  can 
be  described  by  a Gaussian  phase  model,  and  that  these 
phase  distortions  are  the  principal  cause  of  image 
degradation.  These  atmospheric  effects  were  simulated  by 
randomizing  the  phase  of  the  pupil  function  of  the 
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incoherent  imaging  system. 

A series  of  arrays  of  Gaussian  random  noise  variables 
with  unit  variance  was  generated  from  a random  process 
which  was  zero  mean,  unit  variance,  and  uncorrelated. 
Correlated  noise  was  produced  by  filtering  these  arrays  by 
the  technique  described  in  Knox  [15],  with  the  correlation 
length  corresponding  to  a turbulence  cell  size  of  11  inches 
across  a 56-inch  telescope.  Fig.  3 shows  one  of  these 
random  number  arrays.  The  point-spread  function  of  the 
digital  correlation  filter  is  shown  in  Fig.  4.  A 
correlated  noise  array  is  illustrated  in  Fig.  5.  Each 
correlated  noise  array  was  used  as  the  phase  9 ( k ) of  a 
complex  pupil  function  whose  magnitude  equaled  the  aperture 
function  A ( x ) of  Fig.  6.  The  pupil  functions  were  inverse 
Fourier  transformed  and  the  squared  modulus  formed  yielding 
a set  of  45  different  point-spread  functions.  One  of  these 
i s shown  in  Fig  . 7 . 

An  independent  set  of  45  point-spread  functions  was 
similarly  generated.  These  were  convolved  with  the  double 
star  image  of  Fig.  2 to  yield  the  set  of  45  blurred  images 
i (x),  one  of  which  is  shown  in  Fig.  R.  These  blurred 
double  star  images  represent  the  sho r t- ex po sur e photographs 
which  are  then  input  to  the  Knox-Thompson  procedure. 

To  examine  the  sensitivity  of  this  technique  to  sensor 
noise,  we  performed  a series  of  simulations  with  different 
levels  of  sensor  noise  present.  The  sensor  noise  model  is 
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Fig.  3 Random  numbers  (to  be 
used  as  the  Phase  of  the 
Instantaneous  OTF) . 
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Fig.  4 Gaussian  Spatial 
Phase  Correlation  Filter 
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Fig.  5 Spatially  Correlated 
Phase  of  Instantaneous  OTF. 
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based  on  a sem 1 - cl  a s s ic al  approach  to  photon  detection.  In 
this  model,  the  number  of  photoelectrons  released  from  a 
small  region  of  area  A centered  at  a point  (x,y)  in  the 
image  plane  is  taken  to  be  a Poisson  random  variable  with 
mean 

X=[n I ( x,y)/hv]TAA. 

In  this  expression,  n is  the  quantum  efficiency  of  the 
photodetector,  h is  Planck's  constant,  v is  the  mean 
optical  frequency,  r is  the  integration  (exposure)  time, 
and  I(x,y)  is  the  image  intensity  at  point  (x,y). 

It  is  well  known  that  both  the  mean  and  the  variance 
of  the  above  Poisson  random  variable  are  equal  to  x.  Thus, 
the  mean  photoelectron  current  produced  by  a detector  at 
point  (x,y)  is  proportional  to  the  image  intensity  at  that 
point,  but  so  are  the  fluctuations  in  that  current.  If  we 
define  the  SNR  of  the  image  at  point  (x,y)  as  the  ratio  of 
the  square  of  the  mean  current  to  the  variance  of  the 
current,  we  find  that  this  SNR  is  just  X.  We  can  vary  x, 
and  therefore  the  SNR,  by  varying  t (the  integration  time) 
or  aA  (the  area  over  which  the  image  is  averaged). 

An  additional  simplification  results  if  we  assume  that 
the  number  of  photoelectrons  is  large.  In  this  case,  the 
photoelectric  current  will  be  approximately  Gaussian  with  a 
mean  equal  to  its  variance.  We  made  this  assumption  in  our 
simulations.  The  SNR  which  we  ascribe  to  a given  simulated 
image  is  simply  the  maximum  of  the  point-by-point  ratios  of 
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squared  mean  to  variance  as  described  above.  This 
definition  is  somewhat  arbitrary,  of  course,  but  this  is 
not  a serious  problem  since  we  are  interested  primarily  in 
relative  performance  as  we  vary  the  noise. 

Results  of  the  Knox-Thompson  restoration  procedure, 
comparing  both  line  integration  and  1 east- squares 
integration  for  various  SNR  as  defined  aoove,  are  shown  in 
Figs.  9a  through  9y . Although  we  would  like  to  quantify 
the  results  in  some  perceptually  significant  manner,  no 
general  perceptual  distance  metric  has  yet  been  defined  for 
images,  much  less  for  the  way  in  which  phase  errors  might 
affect  image  formation.  We  default  then  to  SSE  measures, 
and  simply  include  illustrations  of  the  various  estimates 
and  associated  reconstructions.  The  WSSE  measure  in  the 
captions  of  the  phase  estimates  refers  to  a weighted  SSE 
measure  where  the  weighting  function  is  the  magnitude  of 
the  object  Fourier  transform.  The  phase  error  measurements 
reflect  the  modulo  2 it  nature  of  the  phase.  For  reference 
the  object  energy  is  2.929  units  and  the  sum  squared  total, 
or  'energy'  of  the  phase,  is  3860.642  units. 

All  the  reconstructions  were  based  on  the  Labeyrie 
magnitude  spectrum  estimate.  To  minimize  artifacts  in  the 
reconstructed  image,  these  magnitude  estimates  were 
digitally  low-pass  filtered  with  a radially  symmetric 
function  which  was  constant  to  25  cycl es/ pi c ture  , and  then 
tapered  to  zero  at  32  cycles/picture  by  means  of  a cubic 
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c)  L i ne- 1 ntegra  1 Object 
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e)  Least-Square  Object 
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d)  Least-Square  Phase 
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9(a-e)  Comparison  of  L i ne- I ntegra  1 and  Least-Square 
Reconstruction  at  5 dB  SNR 


I'v/S** 


foot  w Hui, 


g)  Line-Integral  Phase 
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Comparison  of  L i ne- I ntegra 1 and  Least-Square 
Reconstruction  at  10  dB  SNR 
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Reconstruction  at  20  dB  SNR 
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9(u-y)  Comparison  of  L i ne- I n teg ra 1 and  Least-Square 
Reconstruction  at  40  dB  SNR 


56 


spl i ne  functi on  . 

At  a SNR  of  5 dB,  neither  method  gives  a good 
indication  of  a second  star.  The  least-squares  method  is 
more  localized,  however,  and  does  not  show  the  presence  of 
artifacts  along  the  image  coordinate  axes.  At  10  dB  the 
1 ea st- squares  restoration  still  shows  only  one  star,  while 
the  line  integral  restoration  hints  at  the  possibility  of 
two  stars.  At  20  dB  both  reconstructions  indicate  clearly 
the  second  star.  The  line  integral  version  has  more 
correctly  estimated  the  relative  amplitudes  of  the  stars, 
but  also  has  many  more  artifacts.  At  30  dB  and  40  dB  both 
methods  have  isolated  the  two  stars.  The  1 east- squa res 
versions  are  more  nearly  correct  in  relative  amplitude  and 
are  less  noisy.  (The  SSE  measures  are  somewhat  confusing 
on  the  phase,  where  they  get  better  as  the  estimate  gets 
worse.  Note  also  the  two  SSE  measures  on  the  reconstructed 
objects.  The  first  refers  to  the  error  over  the  entire 
field  of  view,  the  second  to  the  error  within  the  region  of 
support  of  the  star.  For  the  line  integral  object 
reconstructions,  both  measures  behave  erratically). 

Both  types  of  restorations  have  regions  of  negative 
light  intensity.  These  nonphysical  artifacts  are  a result 
of  mi sestimati on  in  the  Fourier  domain,  and  are  quite 
common  in  deconvolution  problems.  They  are  called 
super-blacks  in  the  literature  [55].  In  every  case  the 
peak  negative  light  intensity  is  greater  in  magnitude  in 
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the  line  integral  restorations  than  in  the  least-squares 
restoration,  and  has  a magnitude  of  10  to  20  percent  of  the 
positive  1 ight- intensi ty  peak.  The  restorations  also  tend 
to  be  smeared,  and  occupy  a larger  area  in  the  field  of 
view  than  the  object.  Less  distracting  is  the  fact  that 
the  restorations  are  complex,  instead  of  real.  However, 
the  magnitude  of  the  imaginary  part  is  generally  about  2 
orders  of  magnitude  less  than  the  real  part,  and  is 
entirely  due  to  the  phase  estimate  not  having  perfectly  odd 
symmetry.  In  an  effort  to  further  improve  the 
reconstructions,  a number  of  modifications  were  tried. 

The  first  modification  was  the  inclusion  of  phase 
difference  information  derived  from  the  second-shift  term 


of  the  ACF 

of 

the 

image  transforms.  The 

necessary 

modi fication 

of 

the 

periodic  1 east- squares 

i nteg  rati  on 

algorithm  is  straightforward . We  experimented  with  an 
integrator  which  weighted  differently  the  phase  differences 
from  the  one-  and  two-  shift  terms,  but  found  that  the  WSSE 
of  the  resulting  estimates  was  a monotonical  ly  increasing 
function  of  the  weight  applied  to  the  second-shift  data. 
This  is  consistent  with  Knox's  analysis  and  recommendation 
that  only  one-shift  data  be  used. 

We  also  experimented  with  the  use  of  spatial  domain 
windows  applied  to  each  speckle  photograph.  This  is  a 
common  technique  in  spectral  estimation  for  segmenting  long 
data  sequences  into  manageably  short  lengths.  The  effect 


58 

is  of  course  to  convolve  the  transform  of  each  speckle 
image  with  the  transform  of  the  window  function  before  the 
averaging.  The  effects  on  the  bias  and  variance  of  the 
subsequent  PSD  estimates  are  found  in  Jenkins'  and  Watts' 
text  [56].  We  wanted  to  improve  the  effective  SNR  by 
restricting  ourselves  to  the  brightest  part  of  each  speckle 
image.  We  also  found  in  working  with  the  real  data  that 
the  Knox-Thompson  estimator  seemed  much  less  accurate  than 
our  simulations  would  indicate.  We  hypothesized  that  it 
might  be  quite  sensitive  to  the  violation  of  the 
isoplanatic  assumption,  so  we  wanted  to  restrict  the 
f 1 el d- of- v 1 ew  to  more  nearly  approximate  i sopl anlc i ty . 

Although  the  effects  of  windowing  on  PSD  estimates 
have  been  examined  by  many  researchers,  the  effects  on  the 
Knox-Thompson  estimator  have  not  been  so  catalogued. 
Therefore,  we  examined  the  effects  of  spatial  windowing  on 
the  phase  difference  estimates  by  the  same  simulation 
techniques  previously  described.  A series  of  blurred  stars 
was  created  and  processed  with  each  one  of  three  different 
spatial  windowing  options  to  produce  three  different 
reconstructions.  In  the  first  case  each  speckle  photograph 
was  multiplied  by  a Fourier  (rectangular)  window  of 
dimension  32*32.  The  second  reconstruction  used  a radially 
symmetric  window  made  by  rotating  a Hamming  window  of 
radius  20  about  its  central  axis.  The  third  reconstruction 
used  no  spatial  window,  and  is  that  of  Figs.  do  and  dp. 
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The  example  illustrated  in  Fig.  10  was  performed  with  a 
nominal  SNR  of  30  dB.  In  all  three  cases  the  same 
magnitude  spectrum  estimate  (that  of  the  unwindowed 
version)  was  used,  so  as  to  isolate  the  effects  on  the 
phase  estimation.  As  before,  the  SSE  measures  are  found  in 
the  figure  captions.  Note  that  the  use  of  a Fourier  window 
degrades  the  results,  although  the  Hamming  reconstruction 
is  actually  improved,  at  least  in  terms  of  SSE.  This  is 
not  unexpected  behavior.  Similar  behavior  was  observed  in 
all  other  test  cases  and  in  processing  the  real  data. 

III. 6 Phase  Refinement 

As  mentioned  in  the  literature  review,  the 
Gerchberg- Sax  ton  algorithm  was  created  to  solve  the  phase 
problem  of  electron  microscopy,  where  one  has  knowledge  of 
the  magnitude  of  both  a function  and  its  Fourier  transform 
but  no  knowledge  of  the  phase  in  either  domain.  A flow 
chart  of  the  algorithm  is  found  in  Fig.  11.  The  algorithm 
is  almost  directly  applicable  to  the  problem  of  phase 
estimation  where  one  has  knowledge  of  only  the  magnitude 
spectrum  and  the  region  of  support  of  a real  function  and 
no  knowledge  of  its  phase  spectrum.  However,  this 
relaxation  of  the  required  a priori  knowledge  introduces  a 
serious  non-uniqueness  problem  in  the  reconstructions, 
which  has  been  demonstrated  empirically  by  Feinup  MO]. 
(He  has  some  experimental  evidence  to  indicate  that  the 
problem  may  not  be  as  severe  for  complicated  two 
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d)  Hamming  Window 
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Fig.  10  Comparison  of  the  Effect  of  Different  Spatial 
Windows  on  the  Least-Square  Reconstructions  at  30  dB  SNR 


dimensional  images,  however).  In  this  work  we  have  sought 
to  side-step  the  possible  problem  of  non-uniqueness  by  the 
use  of  the  least-square  Knox-Thompson  estimate  as  an 
initial  approximation  to  the  correct  phase  spectrum.  The 
hope  is  that  it  is  'close  enough'  in  some  sense  to  the  true 
phase  spectrum  to  allow  the  algorithm  to  converge  to  the 
true  phase. 

The  new  algorithm  is  presented  in  flow  chart  form  in 
Fig.  12.  The  only  changes  are  that  instead  of  correcting 
the  magnitude  in  the  spatial  domain,  one  imposes  the 
conditions  of  realness,  non- neg at i v i ty , and  a bounded 
region  of  support  on  the  reconstruction.  To  prove  that  the 
correction  energy  must  decrease  only  a slight  variation 
need  be  made  to  the  geometric  proof  originally  offered  by 
Gerchberg  and  Saxton.  Consider  the  operation  of  the 
algorithm  on  two  discrete  points,  one  in  the  spatial  domain 
and  the  other  in  the  frequency  domain,  as  in  Fig.  13. 
Define  the  squared  error  e as  the  sum  squared  difference 
between  the  original  magnitudes  at  each  discrete  frequency 
and  the  magnitudes  calculated  by  the  algorithm  at 
corresponding  points.  When  e is  zero,  the  algorithm  has 
converged.  Now,  in  the  frequency  domain,  let  the  reference 
magnitude  at  our  arbitrary  point  be  represented  by  the 
measured  magnitude  f Upon  inverse  Fourier 
transformation,  a vector  f-j  is  produced  at  the  arbitrary 
point  in  the  spatial  domain.  If  fi  is  in  the  region  of 
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support  and  if  its  real  part  is  non-negative  the  real  part 
is  preserved.  Otherwise  fj  is  set  to  zero  by  the  addition 
of  a colinear  vector  of  equal  magnitude  and  opposite 
direction.  All  points  in  the  spatial  domain  are  similarly 
corrected  and  the  Fourier  transform  of  the  new  object  is 
taken.  At  our  arbitrary  point  in  the  frequency  domain  the 


new  vector  g^  is  the  sum  of  f^  and  a new  vector  c^  produced 


at  that  frequency  by  the  Fourier  transform  of  the  vectors 


c ^ at  each  point  in  the  spatial  domain.  If  necessary  the 


magnitude  of  g^  is  corrected  to  that  of  f by  the  addition 


of  the  vector  d^,  which  is  colinear  with  g^ . The  algorithm 


now  iterates  until  convergence. 
Parseval's  theorem  that 
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It  is  also  clear  from  examining  the  phasor  plots  of  Fig.  13 
that 

|d2lij:|c2l?o  51>J  <53> 

with  equality  iff  C2  is  colinear  with  f^  . Thus 

(54) 
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with  equality  iff  the  correction  vector  is  colinear  at 


every  frequency.  To  this  point  we  have  shown  that  the 
correction  energy  applied  in  the  frequency  domain  cannot  be 
larger  than  the  correction  energy  applied  in  the  spatial 
domain  just  previously.  Now  , in  the  second  iteration,  we 


have  in  the  spatial  domain  the  vector  f^ ' , the  sum  of 


vectors  g^  and  d^  . It  is  geometrically  obvious  that 
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Thus 


I [e  |f  ; l Id  if— |d  |f 
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(56) 


and  the  correction  signal  energy  in  the  spatial  domain  is 
also  less  than  or  equal  to  the  frequency  domain  correction 
energy  applied  just  prior  to  inverse  Fourier 
transformation.  Hence  the  squared  error  must  decrease  or 
remain  constant  with  each  iteration.  For  the  squared  error 
to  be  bounded  away  from  zero,  it  is  required  that  the 
correction  signal  in  the  frequency  domain  be  colinear  with 
the  estimated  frequency  spectrum  while  being  entirely 
imaginary  in  the  spatial  domain.  This  condition  cannot  be 
fulfilled  because  the  algorithm  forces  the  spatial  estimate 
to  be  real,  and  hence  to  have  an  even  magnitude  spectrum 
and  an  odd  phase  spectrum.  To  be  colinear  the  correction 
signal  must  also  have  odd  phase  symmetry,  and  hence  cannot 
be  Imaginary  in  the  spatial  domain.  (Note  that  this  result 
is  stronger  than  that  presented  for  the  original  electron 
microscopy  problem,  where  certain  phase  symmetries  in  the 
spatial  and  frequency  domains  may  bound  the  correction 
energy  away  from  zero.)  This  result  does  not  necessarily 
imply  convergence  of  the  correction  energy  to  zero,  nor 
does  it  imply  solution  uniqueness,  but  in  our  experience 
this  proviso  is  of  theoretical  rather  than  practical 
interest . 

One  might  wonder  whether,  given  a phase  estimate,  one 
could  effect  a magnitude  solution  by  the  same  algorithm. 


I 


66 


The  answer  is  no,  because  it  can  be  shown  that  the 
correction  signal  energy  cannot  be  guaranteed  to  decrease. 
Imagine  again  the  same  phasor  diagrams,  with  the  same 
initial  signal  estimates,  as  represented  in  Fig.  14.  At 
the  first  return  to  the  frequency  domain,  one  corrects  the 
magnitude  of  g£  by  adding  d£.  Since  62  must  now  be 
perpendicular  to  the  perpendicular  bisector  OP,  it  is 
geometrically  obvious  in  this  case  that  62  > c2  • If 
this  condition  is  satisfied  at  enough  points,  the 
correction  signal  energy  will  increase,  and  the  algorithm 
may  never  converge.  In  fact,  this  behavior  has  been 
observed  in  empirical  tests. 

To  test  the  phase  refinement  technique,  the 
least-squares  phase  estimates  described  in  the  previous 
section  were  used  as  the  initial  approximations  to  the 
phase  spectrum.  In  Figs.  15a  through  15j,  the  refined 
phase  estimates  and  the  resulting  reconstructions  are  shown 
for  each  of  the  various  SNR.  In  each  case,  the  algorithm 
was  stopped  when  the  weighted  correction  signal  energy  was 
reduced  to  .16.  Originally  it  was  intended  to  have  this  be 
.1,  but  the  procedure  was  unable  to  reduce  the  correction 
energy  much  below  .16  for  the  two  lowest  SNR  cases.  In  all 
cases,  the  exact  extent  of  the  object,  or  region  of 
support,  was  assumed  known. 

Phase  refinement  actually  degraded  the  results  for  the 
two  poorest  SNR  cases.  Slight  improvement  for  the  20  dB 
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case  is  apparent,  but  for  the  30  dB  and  40  dB  cases  the 
improvement  is  substantial.  The  refinement  required  29, 
29,  7,  12  and  22  iterations  respectively.  Note  that  as  the 
SNR  gets  poorer,  the  phase  tends  to  be  underestimated.  The 
obvious  conclusion  is  that  the  technique  is  not  a cure-all, 
but  it  can  provide  real  improvement  if  the  magnitude 
spectrum  is  fairly  accurate.  One  would  expect  good  success 
if,  instead  of  the  entire  magnitude  spectrum,  only  the 
frequencies  with  locally  good  SNR  were  used  to  determine 
the  phase  refinement.  This  is  correct,  and  it  forms  an 
important  part  of  the  super- resol ution  procedure  that  is 
discussed  in  the  next  chapter. 
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CHAPTER  IV 

SUPER-RESOLUTION 

IV. 1 Introduction 

One  is  often  confronted  in  a signal  processing 
environment  with  a signal  that  has  been  lowpass  filtered 
with  a resulting  loss  of  information.  One  would  like  to  be 
able  to  recover  this  information  and  restore  the  signal. 
There  is  a strong  basis  in  theory  to  support  the 
plausibility  of  such  recovery  foi  a certain  class  of 
signals,  those  which  have  an  analytic  frequency  spectrum. 
The  difficulty  in  actually  realizing  such  recovery  in  the 
practical  case  has  led  at  least  one  pa i r of  authors  to 
reject  the  concept  as  a 'myth'  [55].  It  is  my  opinion  that 
this  judgement  is  premature  and  that  a substantial  amount 
of  super-resolution  may  be  had  even  in  the  case  of  real 
world  blurs.  This  chapter  discusses  a new  algorithm  for 
achieving  super- resol uti on . It  is  similar  to  the  algorithm 
originally  proposed  by  Gerchberg,  but  in  our  experience  has 
demonstrated  better  convergence  properties  and  much  better 
conditioning  in  the  presence  of  noise.  In  this  chapter  we 
will  present  a review  of  the  Gerchberg  algorithm,  a 
heuristic  analysis  of  the  new  algorithm,  and  the  results  of 
some  simulations  which  compare  the  algorithms  on  the 
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synthetic  blurred  star  data  base. 

IV. 2 Background  and  Analysis  of  the  Algorithm 

The  theory  of  analytic  functions  of  a single  complex 
variable  is  usually  invoked  to  justify  the  attempts  at,  as 
well  as  understand  the  methods  of,  super- resol uti on . In  an 
astronomical  context,  one  is  at  best  limited  to 
d i f f r ac t i on- 1 i m i ted  measurements.  The  OTF  of  a 
d i f f r ac t i on- 1 i m i ted  imaging  system  has  an  absolute  spatial 
frequency  cutoff.  Because  a function  of  finite  support  has 
an  everywhere  analytic  (or  entire)  z transform  [42],  and 
because  an  entire  function  must  be  infinite  in  extent,  we 
can  say  that  d i f f r ac 1 1 on- 1 i m i ted  imaging  always  removes 
spatial  frequency  information,  even  in  the  noiseless  case. 
In  reality,  however,  the  signal  energy  may  be  sufficiently 
concentrated  within  a finite  bandwidth  to  permit  object 
resolution  of  arbitrary  accuracy  with  a finite  spectrum 
[44].  When  this  is  the  case,  the  object  is  said  to  be  well 
resolved. 

When  the  low-pass  filtering  operation  of  the  OTF 
removes  significant  energy  we  say  the  object  is  poorly 
resolved.  The  classical  definition  of  resolution  specifies 
the  minimum  angular  separation  of  two  point  sources  which 
allows  them  to  be  seen  as  distinct.,  for  a given  optical 
instrument  of  circular  pupil  diameter  0.  This  minimum 
separation  was  found  to  be  1 . 22 X/D , where  lambda  is  the 
wavelength  of  the  observed  light..  Oi  Francia  [51]  pointed 


out  the  practical  rather  than  theoretical  nature  of  this 
limit,  however,  and  showed  that  an  infinite  number  of 
object  intensity  distributions  could  produce  an  identical 
image.  By  this  he  demonstrated  the  implicit  use  of  a 
priori  knowledge  on  the  part  of  the  observer  in  determining 
the  resolution  of  an  optical  instrument. 

Harris  [42]  recognized  that  no  two  distinct  objects  of 
finite  angular  size  can  produce  identical  images.  This  is 
because  a "function  of  a complex  variable  is  determined 
throughout  the  entire  z plane  from  a knowledge  of  its 
properties  within  an  arbitrarily  small  region  of 
analyticity"  . Thus  an  entire  function  is  everywhere 
determined  by  its  behavior  over  an  arbitrarily  small 
region,  so  that  an  accurate  measurement  of  a finite  portion 
of  the  spectrum  is  sufficient  to  uniquely  determine  the 
complete  spectrum,  and  hence  the  object,  given  that  it  has 
finite  support.  He  concludes  that  " di f f r ac t i on . . . i mpo se s a 
resolution  limit  which  is  determined  by  the  noise  of  the 
system,  rather  than  some  absolute  criterion." 

While  the  above  reasoning  no  doubt  gives  insight  into 
the  problem,  there  remain  some  unresolved  issues  when  one 
is  constrained  to  computer  solutions.  These  seem  to  me  to 
include  the  following  questions.  First,  is  the  requirement 
that  complete  knowledge  of  the  function  over  some  region  of 
analyticity  satisfied  by  knowledge  of  the  Fourier  transform 
of  the  function,  which  is  at  best  knowledge  along  a line  or 
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contour  in  the  complex  domain?  Second,  is  the  requirement 
satisfied  when  one  has  only  discrete  samples  of  the  Fourier 
transform?  Third,  does  the  proof  of  the  Paley-keiner 
theorem  extend  to  functions  of  more  than  a single  complex 
variable?  In  spite  of  these  questions,  the 
semi-quantitative  analysis  given  here  will  depend  on  the 
properties  of  analytic  functions.  An  analysis  along  the 
lines  of  that  given  by  Youla  [48]  would  be  independent  of 
the  questions  just  raised,  however,  and  would  be  highly 
desi rabl e . 

Since  our  new  algorithm  is  an  adaptation  of  the  one 
due  to  Gerchberg,  a review  of  the  latter  will  help  in 
understand! ng  our  work  and  in  placing  it  in  proper 
perspective.  The  Gerchberg  algorithm  iterates  between  the 
spatial  and  frequency  domains.  In  the  spatial  domain,  the 
reconstruction  is  set  to  zero  outside  the  known  region  of 
support  of  the  object,  and  in  the  frequency  domain  the 
newly  estimated  spectrum  is  always  corrected  to  agree  with 
the  known  portion  of  the  spectrum.  A flow  chart  of  the 
algorithm  is  found  in  Fig.  16.  The  algorithm  is  shown  to 
converge  uniquely  for  noiseless  data.  Gerchberg  recommends 
a procedure  to  limit  the  buildup  of  noise  in  the 
reconstruction  due  to  noise  in  the  'known',  or  measured, 
portion  of  the  spectrum.  This  operation  is  necessary 
because  correcting  the  known  portion  of  the  spectrum  adds 
back  all  the  noise  energy  originally  present,  so  that  the 
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error  energy  due  to  this  noise  must  increase  or  remain 
constant  at  every  iteration.  At  some  point  the  error  due 
to  the  increase  in  this  noise  energy  becomes  greater  in 
magnitude  than  the  decrease  in  error  energy  resulting  from 
refining  the  extrapolated  spectrum.  At  this  point  the 
algorithm  should  be  terminated.  This  result  has  been  shown 
[48]  on  theoretical  grounds  for  the  continuous  case. 
Unfortunately,  the  optimal  number  of  iterations  can  not  be 
known  a priori. 

The  noise  limiting  procedure  is  based  on  the 
observation  that  energy  in  the  converged  reconstruction 
which  lies  outside  the  known  object  extent  must  be  due  to 
the  original  distortion  energy.  Therefore,  the  object  is 
limited  to  its  correct  spatial  extent,  the  transform  taken, 
and  the  new  lowpass  spectrum  is  used  to  start  the  iteration 
all  over  again.  The  procedure  thus  eliminates  that  portion 
of  the  original  spectral  noise  energy  which  lies  outside 
the  known  object  extent. 

The  new  algorithm  corrects  this  deficiency  by  not 
correcting  the  magnitude  spectrum  to  the  measured  magnitude 
spectrum,  so  that  measurement  error  is  not  increased. 
Instead  the  magnitude  and  phase  spectra  are  alternately 
extended  through  the  frequency  space.  The  phase  spectrum 
extension  is  allowed  at  each  step  to  converge  in  a weighted 
SSE  sense,  where  the  weighting  function  is  the  estimated 
magnitude  spectrum.  A flowchart  is  illustrated  in  Fig.  17. 
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As  shown  in  Chap.  III. 6,  no  magnitude  iterations  are 
performed  since  this  is  an  unstable  procedure.  As  each 
frequency  extension  is  made  the  total  energy  in  the 
spectrum  below  the  old  frequency  cutoff  changes.  If  the 
new  magnitude  spectrum  is  simply  scaled  to  agree  with  the 
old  spectrum  (in  whatever  sense  is  appropriate)  below  the 
old  cutoff  frequency,  the  effective  convergence  of  the 
algorithm  is  easily  determined  by  the  amplitude  of  this 
scaling  constant.  In  my  work,  the  algorithm  was  always 
terminated  when  something  in  excess  of  99%  of  the  total 
spectral  energy  was  estimated.  The  phase  estimates 
exhibited  erratic  convergence  behavior  when  frequency 
extrapolation  was  attempted  much  beyond  that  point. 
Although  their  convergence  could  often  be  stabilized  by  a 
smaller  cutoff  frequency  extension,  the  computational  time 
simply  became  excessive. 

There  seem  to  be  two  sources  of  error  in  the 
extrapolated  spectrum  produced  by  the  new  algorithm.  They 
are  induced  by  the  way  in  which  the  algorithm  treats  the 
two  kinds  of  errors  in  the  original  measured  spectrum. 
Following  Gerchberg,  we  call  these  two  kinds  error  energy 
and  distortion  energy. 

We  can  consider  the  initial  lowpass  spectrum  estimate 
to  be  the  sum  of  the  correct  spectrum  and  an  error  spectrum 
which  exists  only  beyond  the  initial  lowpass  cutoff 
frequency.  The  error  spectrum  is  equal  in  magnitude  and 
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Fig.  17.  Flow  chart  of  the  new  nonlinear  super- reso 1 u ti on 
algorithm. 
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opposite  in  sign  to  the  true  spectrum.  Inside  the  cutoff 
frequency,  the  known  spectrum  is  represented  as  the  sum  of 
the  true  spectrum  and  a distortion  spectrum.  Consider 
first  the  case  where  the  distortion  energy  is  zero.  Since 
the  error  spectrum  is  b and- 1 i m i ted  , it  follows  that  its 
inverse  Fourier  transform,  the  error  object,  cannot  be 
spatially  limited.  Thus  some  error  energy  will  be 
preserved  when  the  spatial  constraints  are  imposed. 
Consequently  the  spectrum  of  the  constrained  object  will 
retain  some  portion  of  the  original  error  energy,  some  of 
which  may  be  found  even  in  the  low  frequencies  below  the 
initial  cutoff  frequency.  Since  the  low  frequencies  of 
this  new  spectral  estimate  are  retained  in  the  magnitude 
spectrum  extension,  the  algorithm  has  introduced  into  the 
low  frequency  spectrum  a portion  of  what  was  originally 
error  energy.  To  minimize  the  growth  of  this  error  we 
should  minimize  the  number  of  times  we  extend  the  spectrum. 
This  implies  either  rapid  spectral  extrapolation  or  that  we 
confine  ourselves  to  signals  that  are  nearly  resolved  so 
that  few  extrapol ati ons  are  required.  Because  most  images 
have  a general  lowpass  characteristic,  one  would  expect  the 
error  energy  injected  into  the  low  frequencies  to  decrease 
as  the  cutoff  frequency  is  increased  and  the  remaining 
error  energy  is  decreased. 

Suppose  now  that  distortion  energy  is  present.  It  is 
probable  that  the  spatial  constraints  will  also  preserve 


some  portion  of  this  energy,  as  well.  The  spectrum  of  the 
newly-constrained  object  may  similarly  have  distortion 
energy  present  in  all  frequencies.  To  the  extent  that 
distortion  energy  is  injected  into  the  high  frequencies, 
Parseval's  theorem  indicates  that  the  total  distortion 
energy  must  be  decreased  by  lowpass  filtering,  even  with  a 
higher  cutoff  frequency.  To  minimize  the  distortion  energy 
in  the  final  result  we  want  to  perform  many  filtering 
operations,  which  implies  slow  spectral  extrapolation. 
Changing  the  rate  of  extrapolation  thus  varies  the  relative 
contribution  of  the  distortion  and  error  energies  to  the 
final  extrapolated  estimate.  Although  it  is  possible  that 
a variable  rate  would  Improve  the  algorithm  performance, 
all  the  experimentation  performed  here  had  a fixed  rate  of 
frequency  cutoff  increase.  The  optimum  rate  of  increase  or 
even  Its  existence  in  the  general  case,  is  not  presently 
known . 

The  heuristic  analysis  just  presented  is  linear,  and 
in  the  sense  that  at  every  step  the  algorithm  is  dealing 
with  disjoint  regions  in  the  frequency  domain,  it  is  a 
linear  algorithm.  However  the  possible  importance  of  the 
nonlinearity  introduced  by  the  separate  extension  of  the 
magnitude  and  phase  is  not  yet  known.  In  summary,  the 
original  Gerchberg  algorithm  is  able  in  the  limit  to 
decrease  the  error  energy  to  zero  at  the  expense  of 
Increasing  the  total  distortion  energy.  The  new  algorithm 
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sacrifices  the  possibility  of  perfect  reconstructions  of 
noise-free  data  in  exchange  for  stable  behavior  of  the 
distortion  energy. 

IV. 3 Super- resol uti on  of  Simulated  Data 

I have  examined  by  simulation  the  behavior  of  the 
algorithm  as  the  rate  of  increase  of  the  cutoff  frequency 
is  changed,  although  the  examination  was  certainly  not 
exhaustive.  The  results  on  the  S dB  SNR  case  are 
summarized  in  Table  III,  where  it  appears  that  the  optimum 
increase  is  about  5 discrete  frequencies  at  each  step,  at 
least  for  this  data.  This  same  rate  of  increase  was  used 
in  other  simulations  where  the  effects  of  SNR  and  initial 
lowpass  cutoff  frequency  were  investigated.  In  all  the 
simulations,  the  lowpass  filter  is  a radially  symmetric 
Fourier  lowpass  filter.  The  phase  convergence  criterion 
was  that 

e>ni0(C)|{  ♦(k)(u)-£(k’1)(u))2} 

u 

Since  our  heuristic  analysis  indicates  that  some  error 
must  always  remain  in  the  solution,  and  may  in  fact 
increase  in  the  'known'  portion  of  the  spectrum,  it  was 
necessary  to  observe  the  seriousness  of  this  characteristic 
in  a practical  setting.  A number  of  reconstructions  based 
on  lowpass  filtering  the  non-noisy  data  at  various  cutoff 
frequencies  were  effected.  They  are  shown  together  with 
the  initial  estimates  and  their  super- resol ved  spectra  in 
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Figs.  18  a)  through  18  t)  . In  every  case,  e=.l.  The 
reconstructions  are  quite  stable  with  decreasing  initial 
frequency  cutoff,  although  it  is  curious  to  note  that  the 
last  estimate  has  a slightly  better  SSE  than  the  third 
estimate.  The  recon struc ti on s required  10,  15,  18,  and  22 
iterations,  respectively.  All  of  the  reconstructions  have 
errors  less  than  2.5%.  It  appears  that  the  error 
introduced  into  the  'known'  portion  of  the  spectrum  is 
relatively  small  and  well  behaved. 

Figs.  19  a)  through  19  y)  show  the  super- resol ution  of 
the  least-squares  phase  Knox-Thompson  reconstructions  of 
chapter  III.  It  was  necessary  to  use  a lower  and  lower 
initial  cutoff  frequency  as  the  SNR  got  poorer,  so  as  to 
maintain  a reasonably  accurate  initial  estimate.  Note  that 
even  in  the  5 dB  SNR  case  an  excellent  reconstruction  is 
produced.  In  most  cases  the  peak  magnitudes,  while  about 
correct  relative  to  each  other,  are  in  error  absolutely. 
The  information  in  the  peaks  resides  almost  entirely  in  the 
highest  frequencies,  where  less  than  one  percent  of  the 
signal  energy  is  concentrated.  Mi sestimati on  in  these 
frequencies  is  thus  easily  detected  in  the  spatial  domain. 
(The  final  estimates  have  been  scaled  to  have  total  energy 
equal  to  that  of  the  true  object).  None  of  the  solutions 
required  more  than  33  iterations. 

The  performance  of  the  procedure  seems  to  depend  not 
only  on  the  amount  of  noise  in  the  initial  estimate  but 
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also  on  its  structure,  so  that  an  apparent  increase  in 
local  SNR  may  not  necessarily  imply  a more  accurate 
reconstruc ti on , The  spectral  estimates  are  smoothed  by  the 
procedure,  as  shown  in  the  illustrations.  Particularly 
striking  is  the  fact  that  the  distortion  energy  in  the 
Initial  magnitude  spectrum  has  been  decreased  in  the 
corresponding  lowpass  portion  of  the  extrapolated  spectrum 
to  between  It.  and  lit  of  its  original  value.  However,  in 
extrapolating  real  frequency  data,  sit ua t ions  might  arise 
where  one  would  prefer  to  replace  the  final  low  frequency 
spectrum  with  the  measured  spectrum  to  avoid  as  much  as 
possible  this  smoothing. 

As  an  Illustration  of  the  superiority  of  the  new 
algorithm  a reconstruction  of  the  5 dB  case  using 
Gerchberg’s  algorithm  is  Included.  The  same  initial 
lowpass  filter  was  used,  and  the  phase  was  refined,  so  that 
for  both  algorithms  the  'known'  portion  of  the  spectrum  was 
Identical.  In  Figs.  20  a-c)  a solution  obtained  by 
allowing  150  iterations,  with  the  noise  limiting  procedure 
on  the  101st  iteration,  is  shown.  This  reflects  the  common 
usage  of  the  algorithm.  In  Figs.  20  d-f)  the  optimal 
number  of  iterations  was  allowed,  that  Is,  29  Iterations 
with  the  noise  limiting  procedure  on  the  15th  iteration. 
(It  is  clear  that  this  determination  of  the  optimum  number 


was  possible  only  because  the  answer  was  alreaoy  known). 
In  both  cases  very  inferior  reconstructions  were  produced. 
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This  disparity  lessens  as  the  quality  of  the  initial  data 
improves,  but  the  superior  convergence  properties  of  the 
new  algorithm  give  it  the  advantage  in  these  cases  as  well. 
As  the  noise  In  the  data  decreases,  the  amplification  of 
the  noise  by  Gerchberg's  algorithm  does  not  become  a 
serious  problem  until  the  error  spectrum  energy  is  reduced 
to  a lower  level.  Thus  the  optimal  number  of  iterations 
increases.  It  is  not  uncommon  to  use  several  hundred 
iterations  to  approximate  convergence.  The  new  algorithm 
seems  always  to  converge  in  about  the  same  number  of 
Iterations  for  a given  initial  cutoff  frequency,  nearly 
independently  of  SNR. 

To  test  the  sensitivity  of  the  reconstructions  to 
inaccuracies  In  determining  the  region  of  support  a number 
of  mask  sizes  was  used.  The  true  region  of  support  was  18 
pixels.  Illustrated  In  Figs.  21  a)  through  21  f)  are 
reconstructions  based  on  initial  non-noisy  data  using  masks 
of  area  24  and  32  pixels,  respectively,  assuming  an  initial 
cutoff  frequency  of  7 cycles  per  picture,  as  in 
Figs.  18  p-q)  . The  resulting  reconstructions  have 
substantial  error,  although  the  binary  star  character  of 
the  object  is  still  clearly  visible.  Of  course,  the 
changes  in  the  lowpass  magnitude  spectrum  are  indicative  of 
the  error  In  estimating  the  region  of  support.  We  have 
also  observed  that  the  phase  convergence  can  be  adversely 
affected  if  the  region  of  support  is  poorly  known.  The 
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Fig.  21  Super-Resolution  of  Noise-Free  Data  Using 
Inaccurate  Spatial  Constraints 
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importance  of  an  accurate  estimate  of  the  region  of  support 
is  certainly  indicated,  although  the  degradation  with 
increasing  mask  size  is  graceful. 

In  super- resol vi ng  real  data,  it  might  be  difficult  to 
accurately  estimate  the  region  of  support  based  on  only  the 
estimated  autocorrelation  of  the  object,  which  will  itself 
have  been  low-pass  filtered.  It  is  quite  possible  that  a 
priori  knowledge  may  be  needed  to  determine  the  object  size 
and  shape.  An  alternative  would  be  simply  to  iteratively 
approximate  the  region  of  support  until  one  determined 
which  region  minimized  the  change  in  the  'known1  portion  of 
the  magnitude  spectrum,  subject  to  whatever  knowledge  is 
available  concerning  the  object  itself. 
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CHAPTER  V 


RESTORING  TURRULE NC E -B LUR RED  STARS 

V. 1 Description  of  Data  Base 

The  techniques  described  in  the  preceding  chapters 
have  been  applied  to  turbul ence- bl ur red  images  of  the  stars 
a-ORI  (Betelgeuse)  and  a-AUR  (Capella).  The  point  star 
reference  used  to  estimate  to  mean-square  OTF  was  B-ORI. 
Betelgeuse  is  a red  super-giant  star  with  the  largest 
apparent  angular  diameter  of  any  single  star  in  the 
heavens.  Capella  is  a binary  star  whose  maximum  interstar 
distance  is  less  than  .06  arc-seconds,  which  is  comparable 
to  diameter  of  Betelgeuse.  Capella  is  particularly 
interesting  in  that  it  is  unresolved  in  our  data. 

Our  data  base  consists  of  180  images  of  each  star. 
The  observations  were  made  with  the  160-inch  Mayall 
telescope  at  the  Kitt  Peak  Observatory  in  Arizona,  by  S.P. 
Worden  and  B.S.  Baxter,  on  the  nights  of  June  13  and  14, 
1976.  The  optical  interference  filters  had  a central 

o o 

wavelength  of  6500  A and  a bandpass  of  20  A.  Exposure 
times  were  about  50  ms,  with  approximately  one  second 
between  exposures  to  allow  the  atmospheric  point  spread 
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processed  on  an  image  intensifies  sampled  and  the  image 
Intensity  was  digitized  to  12  bits  on  a 256*256  grid.  The 
scanning  aperture  was  square  with  width  100  microns  and  the 
sampling  Interval  was  50  microns.  The  resulting  scale  is 
.01  arc- seconds/grid  Interval,  for  a square  field  of  view 
2.56  arc-seconds  on  a side.  Since  at  this  wavelength  the 
absolute  frequency  cutoff  of  the  telescope  is  about  27 
eye  1 es/ arc  second , or  about  69  eye  1 es/ pi c ture , the  speckle 
images  are  oversampled  by  a factor  of  about  1.85. 

V.2  Betelgeuse  reconstruc ti on 

A typical  speckle  image,  this  one  of  Betelgeuse,  is 
shown  In  Fig.  22.  The  corresponding  Labeyrie  magnitude 
estimate  is  shown  in  Fig.  23,  while  the  1 east- squares  phase 
estimate  derived  from  processing  45  unwindowed  speckle 
Images  is  shown  in  Fig.  24.  The  corresponding  spatial 
image  is  in  Fig.  25.  (This  latter  has  been  Interpolated  by 
surrounding  the  spectrum  with  data  points  equal  to  zero. 
The  effect  is  to  increase  the  resolution  of  the  DFT,  as  is 
well  known.  In  this  case,  the  zero-padded  array  was  16 
times  larger  in  area.  Increasing  the  DFT  resolution  by  a 
factor  of  4.)  The  failure  of  the  technique  is  evident. 

The  reconstruction  of  Figs.  26  a-c)  is  based  on  the 
same  speckle  images,  this  time  using  a spatial  window.  The 
window  was  a rotated  version  of  a Hamming  window  of  radius 
40  points.  The  bulk  of  the  degradation  is  now  removed, 
although  the  object  is  insufficiently  resolved.  Much  of 
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Fig.  26  Betel geuse  Reconstructions  Obtained  from  Hamming 
Windowed  Speckle  Photographs 


the  problem  Is  due  to  the  rather  poor  magnitude 
characteristics  of  the  final  digital  lowpass  filter  which 
Is  applied  to  eliminate  the  noise  present  In  frequencies 
above  the  diffraction  limit.  To  minimize  the  annoying 
artifacts  of  ringing  Induced  In  the  image  by  the  filtering 
operation,  one  Is  constrained  to  filters  with  gentle 
roll-off  characteristics.  Hence,  additional  smoothing, 
beyond  that  produced  by  the  telescope,  is  present. 

The  results  of  Initial  phase  refinement  at  frequencies 
less  than  25  cycl es/ arc- second , followed  by 
super- resol utl on , are  shown  In  Figs.  27  a-d),  and  in 
Fig.  28.  The  need  for  a final  digital  lowpass  filter  has 
been  eliminated,  since  the  algorithm  forces  the  solution  to 
be  free  of  ringing.  This  Is  perhaps  the  primary  effect  of 
super- resol utl on  on  this  reconstruction,  since  estimating 
the  highest  frequencies  Increased  the  total  energy  by  only 
about  3%.  The  Image  shows  little  evidence  of  fine  surface 
structure,  In  contrast  to  the  result  of  McDonnell  and  Bates 
[7], 

As  discussed  In  chapter  IV,  there  are  some  Internal 
evidences  that  can  be  examined  to  evaluate  the  accuracy  of 


the  reconstruction. 
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frequencies  within  the  passband  of  the  telescope.  This  is 
consistent  with  the  results  obtained  In  the  simulations 
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using  noise-freo  data.  Second,  the  super- resol ved 
magnitude  spectrum  decays  nicely  with  increasing  frequency, 
and  shows  no  suspect  structure  at  the  high  frequencies. 
Third,  the  diameter  of  the  reconstruction,  which  is 
.065  + .005  arc-seconds,  agrees  very  well  with  the 
autocorrel ati on  estimate  derived  from  the  unwindowed 
speckle  images,  which  has  a diameter  of  .116  arc-seconds. 
The  reconstruction  has  been  examined  by  astronomers,  who 
indicate  that  the  average  radial  intensity  profile  of 
Fig.  27  d)  is  In  excellent  agreement  with  astrophysical 
models  of  the  star.  Fourth,  the  peak  negative  light 
Intensity  in  the  Interpolated  image  Is  only  about  IT  of  the 
peak  positive  light  intensity.  Finally,  because  the 
spatial  mask  was  square  of  width  7 pixels.  It  cannot  be 
claimed  that  the  spatial  mask  produced  the  symmetry  of  the 
resulting  reconstruc ti on . The  fact  that  the  reconstruction 
has  circular  symmetry  when  the  spatial  mask  approximated  so 
poorly  the  region  of  support  is  due  undoubtedly  to  the  fact 
that  the  spectrum  was  very  nearly  resolved  by  the 
telescope.  Thus  the  primary  effect  of  the  processing  was 
to  provide  an  accurate  phase  estimate. 

V.3  Capella  reconstruction 

The  reconstruction  of  Capella  was  less  direct  than 
that  of  Betelgeuse  because  of  the  difficulty  in  determining 
the  size  and  spacing  of  the  spatial  masks.  The  Labeyrie 
autocorrelation  estimate  is  so  severely  filtered  that  only 
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the  relative  orientation  of  the  two  stars  and  a very  loose 
bound  on  the  area  of  the  region  of  support  could  be 
obtained  from  it.  The  sizes  and  spacing  of  the  final  mask 
were  arrived  at  by  trial  and  error  and  represent  the  best 
fit  to  the  measured  frequency  spectrum.  I used  as  a 
goodness-of- f i t criteria  the  following  parameters:  1)  the 
difference  in  energy  between  the  initial  spectrum  and  the 
corresponding  lowpass  portion  of  the  extrapolated  spectrum, 
2)  the  location  and  magnitude  of  the  first  spectral  peak 
and  valley,  and  3)  the  negativity  in  the  reconstructed 
image  after  it  had  been  interpolated  by  a factor  of  4.  It 
was  necessary  in  the  processing  to  obtain  finer  resolution 
of  the  image  In  order  to  adequately  approximate  the  size  of 
the  star  discs.  The  measured  frequency  spectrum  was  low 
pass  filtered  with  a separable  Fourier  window  with  cutoff 
frequency  of  27  cycles/arc-second.  The  final  spatial  mask 
used  two  radially  symmetric  masks  of  diameter  5 pixels 
spaced  diagonally  from  each  other  8 pixels  over  and  up. 
This  corresponds  to  an  interstar  distance  of  .028 
arc-seconds,  and  a common  star  diameter  of  about  .0125 
arc- seconds . 

The  spectra  illustrated  in  Figs.  29  a-b)  were  obtained 
from  45  Hamming  windowed  speckle  photographs.  The  profile 
of  the  magnitude  spectrum  taken  along  the  antidiagonal 
makes  it  easier  to  spot  the  telescope  cutoff  frequency,  at 
about  69  eye  1 es/ pic ture  (27  cycl es/arc- second) , as  shown  in 


!■  j. 

] I 

\ 1 


* 


1 


Fig.  29  c).  Note  that  the  phase  estimate  has  failed  except 
at  the  lowest  frequencies.  The  Labeyrle  autocorrelation 
estimate  Is  shown  In  Fig.  29  d) . The  resolution  Is  so  poor 
that  a meaningful  image  cannot  be  formed  from  this  data. 
The  super- resol ved  spectrum  is  shown  in  Figs.  30  a-c).  The 
energy  in  the  spectrum  is  now  4.18  times  the  originally 
measured  energy,  and  the  highest  extrapolated  frequencies 
are  now  at  147  eye  1 es/ arc- second  (resolution  equivalent  to 
an  800-inch  telescope).  The  first  spectral  valley  and  peak 
are  in  essentially  the  same  locations,  and  the  SSE  between 
the  measured  spectrum  and  the  corresponding  portion  of  the 
extrapolated  spectrum  is  2.7%.  The  smoothing  effects  of 
the  algorithm  are  especially  evident  In  the  antidiagonal 
profile  of  the  extrapolated  spectrum.  The  reconstructed 
object.  Interpolated  by  a factor  of  16,  Is  shown  as  a 
line-drawing  in  Fig.  30  d)  and  as  an  Intensity  object  in 
Fig.  31.  The  object  had  a maximum  negative  value  2.85%  of 
the  peak  positive  Intensity.  The  two  stars  have  an 
Intensity  ratio  of  1 : -638.  This  value  is  given  In  [60]  as 
1 : .632  . The  reconstruction  thus  appears  to  be  reasonably 
accurate  . 


It  Is  obvious  In  inspecting  the  reconstruction  that 
the  spatial  masks  have  Interfered  with  estimating  the  star 
discs.  Although  the  dim  star  diameter  is  about  correct, 
the  bright  star  diameter  has  definitely  been 
underestimated.  In  [61]  it  is  reported  that  the  bright  and 


108 


dim  stars  have  diameters  of  13  and  7 times  the  solar 
diameter.  It  thus  appears  that  even  this  artifact  Is 
consistent  with  current  astrophyslcal  measurements. 
Because  the  spatial  mask  has  In  this  case  such  a strong 
Influence  on  th'e  reconstruction,  one  would  like  to 
exhaustively  explore  the  effects  of  variations,  such  as 
using  masks  of  slightly  different  sizes  for  the  two  star 
discs.  It  Is  likely  that  better  results  can  be  obtained  by 
such  experiments.  Unfortunately,  really  fine  spatial 
resolution  would  require  huge  arrays,  and  is  entirely 
beyond  the  practical  capabilities  of  our  computing 
fac llltles. 
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a)  Labeyrie  Magnitude 
Spectrum  Estimate 


b)  Least-Squares  Phase 
Estimate 


c)  Ant i -Di agona  1 Profile  of 
Magnitude  Spectrum  to 
70  cycles/picture 


d)  Labeyrie  Auto-Correlation 
Estimate 


29  Capella  Reconstructions  Obtained  from  Hamming 
Windowed  Speckle  Photographs 
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c)  Anti -di agona  1 Profile  of 
Magnitude  Spectrum  to 
70  cycles/picture 


d)  Intensi  ty  Object 


Fig.  30  Super- Resol ved  Capella  Reconstruction 


CHAPTER  VI 

SUMMARY  AND  CONCLUSIONS 

As  claimed  in  the  introduction,  we  have  presented  a 
method  of  obtaining  reliable  phase  spectrum  estimates  for 
objects  blurred  by  atmospheric  turbulence.  The  method  is 
based  on  the  Knox- Thompso n phase  difference  estimator 
which,  without  modification,  seems  unreliable  and  overly 
sensitive  to  noise.  We  have  also  presented  a new 
super- resol uti on  algorithm  which  appears  to  have  the  best 
stability  and  convergence  properties  of  any  such  algorithm 
proposed  to  date.  Both  methods  have  been  applied  to 
simulated  and  real  blurred  stars  with  substantial  success. 

As  with  most  other  new  techniques,  however,  a number 
of  issues,  both  pragmatic  and  theoretical,  need  further 
study.  The  first  obvious  question  relates  to  the  accuracy 
of  the  reconstructions  of  Betelgeuse  and  Capella.  There  is 
some  recent  evidence  [52]  which  suggests  that  our  speckle 
photographs  are  not  of  the  highest  quality.  Specifically, 
one  would  expect  better  results  if  exposure  times  were 
about  10  ms,  rather  than  50  ms.  This  would  necessitate  the 
use  of  wider  bandpass  optical  interference  filters  to 

o 

permit  adequate  exposure  of  the  speckle  images,  say  100  A 

o 

instead  of  the  20  A of  our  data.  Certainly,  a study  of  the 
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effects  of  the  various  exposure  time-filter  bandwidth 
tradeoffs  should  be  made  for  astronomical  applications. 
Furthermore,  solid-state  cameras  are  now  available  which 
record  a sampled  and  digitized  image  directly  onto  magnetic 
tape.  This  promises  important  reductions  in  sensor  noise 
level  in  the  data.  Consequently,  we  feel  that  although  our 
reconstructions  are  probably  the  best  to  date,  they  cannot 
be  considered  as  the  last  word  on  the  a s trophy s i c al 
structure  of  these  stars.  Rather,  they  indicate  the  power 
of  these  new  processing  techniques.  A possible  empirical 
check  on  these  techniques  would  be  to  reconstruct  an  object 
for  which  we  have  a satellite  photograph  uncorrupted  by 
atmospheric  turbulence.  The  Galilean  moons  of  Jupiter  are 
possible  candidates.  If  our  results  are  confirmed  by 
further  research,  an  exciting  astronomical  application  of 
the  new  super- resol ution  algorithm  would  be  to  data  from 
the  orbital  telescope  now  under  construction. 

Another  research  area  is  the  application  of  these 
methods  to  spectral  estimation,  bandwidth  compression,  and 
blind  deconvolution  problems  for  noth  one-  and 
two-dimensional  signals.  In  this  latter  problem,  one 
requires  a prototype  spectrum  so  as  to  be  able  to  estimate 
the  spectrum  of  the  blur,  which  is  then  removed  from  the 
blurred  signal  by  some  inverse  filtering  technique.  We 
remind  the  reader  of  the  importance  of  an  'adequate* 
initial  phase  spectrum  estimate,  which  may  inhibit  the 
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application  of  the  supe r- reso 1 u t i on  technique  to  certain 
probl em  s . 

There  is  no  doubt  that  our  present  theoretical 
understanding  of  these  techniques  is  incomplete.  For 
example,  rather  detailed  studies  of  the  bias  and  variance 
characteristics,  and  how  these  are  modified  by  analysis 
windows,  have  been  made  for  PSD  estimators.  For 
completeness,  similar  studies  should  be  made  for  the 
Knox-Thompson  estimator.  However,  I feel  an  even  more 
important  development  would  be  the  creation  of  a general 
analytical  framework  in  which  to  cast  the  supe  r- reso  1 ut.  i on 
algorithm  (which  includes  the  phase  refinement  algorithm). 
Youla's  work  T48]  is  probably  the  foundation  of  such  a 
framework.  Our  algorithm  appears  to  be  a special  case  of 
the  regularizors  discussed  in  that  paper.  The  analysis  is 
complicated  by  the  non-linearities  of  the  algorithm, 
however.  Hopefully,  such  an  analysis  would  provide  a 
quantitative  description  of  the  effects  of  noise  and  the 
rate  of  frequency  ex tr a po 1 a 1 1 on  on  the  final  estimates.  In 
fact,  such  an  analysis  may  indicate  that  the  algorithm 
presented  here  is  not  optimal.  This  would  indicate  the 
possibility  of  even  better  results  than  those  presented 
here.  As  a practical  matter,  however,  the  simplicity  of 
the  implementation  and  the  rapidity  of  convergence  of  this 
new  algorithm  may  make  it  a processing  tool  of  importance. 
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